Lidar Modality Handbook

This application area is currently ranked as mature. The discussion below lists additional features and capabilities that can still be addressed.

The existence of these limitation is primarily a result of funded research priorities. Our past and current research has not placed a priority on resolving this current set of limitations. As researchers, we are always seeking opportunities to to address known limitations if funding is available.

To maximize the flexibility in modeling end-to-end LIDAR systems, a full simulation is divided into two, discrete stages. Stage 1 is the "Radiometric Simulation" that (for each pulse) computes the temporal flux received by each pixel during a user-specified listening window. That result is then input into Stage 2, or the "Detector Model and Product Generation" stage, where the output of Stage 1 is coupled with a separate, back-end detection model to convert the temporal photon arrival times into common data products (e.g. geo-located point clouds).

The "Detector Modeling and Product Generation" tool models (a) the unique characteristics of the detector type employed in the receiver to generate time-of-flight measurements and (b) converts those time-of-flight measurements into 3D points within a 3D point cloud (see the figure below). The DIRSIG model is distributed with multi-return Linear-Mode (Lm), single-return Geiger-Mode (Gm) and multi-return photomultiplier tube (PMT) detector models that will be discussed in the later sections of this document.

This two-stage approach provides a lot of flexibility for doing detector performance studies because the detection and product generation can be repeatedly performed very rapidly once the intermediate radiometric data product has computed. The details of these two stages are discussed in later sections of this document.

At the heart of a DIRSIG LIDAR simulation is a flexible and detailed description of the acquisition platform. Like the passive systems modeled with DIRSIG, an active system can be placed at a fixed location (for example, a tripod) or on a moving platform (for example, a ground vehicle, aircraft or spacecraft). The details of platform positioning are beyond the scope of this document. Instead this document will attempt to describe the detailed transmitter and receiver characteristics available to the user. These include (but are not limited to):

The addition of an active, laser radar capability to the DIRSIG4 model was accomplished by the addition of a suite of new objects to the existing radiometry framework. In general, the model is designed to predict the returned fluxes from the scene as a function of time with respect to the shooting of the source laser. The specific challenges of this imaging model were driven by the requirement to predict the received photon counts as a function of space and time. The photon flux arriving at a LIDAR system often approaches discrete photon events due to the low amount of backscattered radiation and can prove difficult for many traditional Monte-Carlo ray tracing techniques. Additionally, the temporal structure of these returns is driven by the spatial structure of the scene and the total travel time of arriving photons accrued during multiple bounce and scattering events within the scene. Analytical, statistical, and existing passive radiometry solvers were found to be insufficient in many instances, particularly for low flux situations. Thus, a new approach was sought.

We hybrid approach we implemented is called "Direct Return and Photon Mapping" or "Direct+PM", for short. This hybrid approach leverages a modeling technique called "photon mapping" (Jensen, 2001). ^[2]. The "direct return" aspect entails backward ray-tracing from the receiver focal plane(s) into the scene. At each intersection, a direct (single bounce) contribution is computed by ray tracing to the transmitter (source). This calculation incorporates surface and/or volume optical properties and the characteristics of the source (spatial/angular intensity, etc.).

The multiple bounce (or multiple scattering) contribution for the location employs a modification to the photon mapping approach. One reason that this technique was selected was because photon mapping has been demonstrated to be applicable to traditional solid geometry reflective illumination and scattering and absorption by participating mediums, particularly in multiple bounce and multiple scattering cases. The photon mapping approach is a hybrid of traditional forward and backward Monte-Carlo ray tracing techniques. In this two-pass method, source photon "bundles" (a set of literal photons that are assumed to travel together) are shot from a source into the scene using forward ray tracing during the first pass and then collected using a backward ray tracing during the second pass. The collection or rendering process utilizes the events recorded during the first pass to calculate the sensor reaching radiance. For the purpose of active laser radar applications, some modifications to the basic photon mapping treatment were made including the tracking of total travel time and a literal photon counting process.

The issue with the photon mapping approach used in DIRSIG4 was that the user was required to correctly setup key parameters for the algorithm to work correcly. Incorrect configuration of the photon mapping parameters could yield incorrect link budgets, noisy returns, etc. One of the overarching goals of DIRSIG5 was to eliminate the dependence on the user correctly configuring the complex algorithms used by the radiometry solvers, photon mapping, etc. Hence, In DIRSIG5, the radiometry core was united under the path tracing numerical radiometry method that provides the user with a single variable to control the radiometric fidelity (how many samples or paths per pixel). DIRSIG5 still incorporates a bi-directional propagation approach but it is fully automated and does not require user configuration. The DIRSIG5 approach models the laser pulse as a geometric mesh that can be propagated, distorted, refracted, etc. The path tracing method can query these pulses and incorporate the energy into the solution.

The contents of a BIN file can be quickly summarized by an analysis tool that is included with DIRSIG. This tool can be accessed on the command-line as the program bin_analyze or via the Detector Model graphical user interface via the Analyze button on the Input tab.

The program will print out the contents of the various headers in the BIN file and will include basic statistics for each pulse. Example output of a BIN file containing a single pulse that was created with a single pixel receiver is shown below:

The bin_analyze tool also has two options that are currently only accessible via the command-line.

Including --output_profiles=true in the arguments will tell the tool to generate an ASCII/Text file containing the average return profile for all the pixels in the array. A separate file will be generated for each pulse in the BIN file using a naming scheme of profile_N.txt, where N is the pulse number starting at 0. The file contains six columns:

The plot below was created from a "profile" for a single pulse for a LIDAR receiver containing a single-pixel array by plotting column #1 vs. column #3. Note that because this was a single-pixel array, column #3 and #4 had the same values.

Including --output_grids=true in the arguments will tell the tool to generate an ASCII/Text file containing the temporally summed photons for each pixel in the array. A separate file will be generated for each pulse in the BIN file using a naming scheme of pulse_N.txt, where N is the pulse number starting at 0. The file contains three columns:

The visualization below was created by plotting an output "grid" file for for a single pulse onto a 128 x 128 pixel array. The data was plotted using the X and Y values (columns #1 and #2) and using the total number of photons (column #3) to drive the color.

The BIN file produced by the "Radiometric Simulation Stage" is input into the "Detector Modeling and Product Generation Stage" by starting the LIDAR Detector Model from the Lidar Tools sub-menu on the Tools menu. The Input tab allows the user to add one or more input BIN files to the model (multiple BIN files might be produced by a single simulation if separate files are created for each task or each capture (pulse)).

A Linear-mode Avalanche Photo-Diode (LmAPD) device is operated just below the device’s breakdown voltage. Under this operating mode, the device is referred to as a "linear" detector because the generated photocurrent is proportional to the received optical power (it should be noted that this relationship isn’t necessarily perfectly linear at all bias voltages). Generally, LmAPD systems have been employed in higher link budget applications, but as fabrication improvements are made the potential for achieving photon counting sensitivity is increasing.

Most Linear-mode (Lm) systems are commercial and as a result specific details about how they work and how they can be modeled isn’t openly published. In practice, the output of an Lm detector is connected to real-time analysis hardware that is "looking" at the signal output as a function of time. The common method used in the real-time analysis is a "constant fraction discriminator" (CFD) rather than a simple threshold. This allows the hardware to detect peaks at various magnitudes.

The linear-mode detection model included with DIRSIG is a basic CFD based detection and timing scheme. The implementation of the CFD in the DIRSIG linear-mode detection model is as follows:

To illustrate this process, a plot of the these various components are shown below (the X axis has been zoomed to emphasize details). The original received waveform is plotted in red, the shifted waveform is plotted in green and the difference waveform is plotted in blue. The resulting triggers (where the blue difference waveform crosses zero with a positive slope) are plotted as the vertical magenta lines. Note how these trigger lines also capture the magnitude at these trigger times.

The LmAPD detector model is selected from the Detector Type menu on the Detector tab in the Lidar Detector Model. This provides the interface for the user to set the various parameters describing the detector’s characteristics.

A Geiger-mode Avalanche Photo-Diode (GmAPD) device is a photon counting device that employs a solid state photo-diode that is biased just above the breakdown voltage of the device. A single photo-electron can trigger a breakdown event, generating an avalanche that can directly trigger a timer. However, the device must then be reset after each breakdown event in order to detect a new arrival event, and this process is slow enough that it cannot be accomplished within a single listening window. Hence, the current generation of GmAPD devices can only make a single range measurement per pulse.

The Geiger-mode (Gm) detection model used is an implementation of the model described by O’Brien and Fouche ^[3]. This approach includes treatment of noise, photoelectron generation and probabilistic detection. The probability of the device firing in bin i is:

\( P(i) = \exp (-fN) [1 -\exp(-S(i) -w)] \)

where f is the fraction of the bins in the range window before the current bin, N is the noise rate per pulse (photoelectron events from passive background flux and dark count across all bins in the listening window), S is the received active signal (in photoelectrons) and w is the noise rate per bin (photoelectron events per bin from passive background flux and dark count). The exp(-fN) term represents the probability that the detector has not been triggered in a previous time bin (remember that a GmAPD device can only trigger once per pulse). The [1-exp(-S(i)-w)] term represents the probability of a trigger in the current time bin due to either signal or noise.

The photon detection efficiency (PDE) and dark count rate (DCR) are incorporated as follows. The signal (S) in photoelectrons is computed from the number arriving photons at the surface of the device multiplied by the PDE. The noise rate per pulse (w) is computed as the dark count rate (DCR) and the background (passive) flux times the length of the listening window (range gate). The noise rate per bin (N) is compute as the noise rate per pulse (w) divided by the number of time bins per listening window.

The equation to compute probability of a trigger in a given time bin can be used to determine at what time the detector will trigger using a Monte-Carlo process. Starting in the first time bin, a random value can be drawn and compared to the probability of that bin. If the random value is lower than the probability, then a trigger occurs in that bin. If not, this process is repeated with the next time bin.

Although this Monte-Carlo approach is correct, it requires a random number to be generated for every time bin. An alternative method is to compute cumulative probabilities across time bins and then draw a single random value to determine during which time bin the detector will fire in. This approach is more attractive computationally because a significantly smaller number of random numbers need to be generated.

The modeling technique cumulative probability approach is as follows:

The current generation Geiger-mode devices can only trigger once while it is armed (during the listening window). Multiple returns with GmAPD system are obtained by imaging the same location multiple times. Because the detection is a Poisson process, a detector illuminated with the exact same waveform will not trigger at the same time on subsequent measurements. To illustrate this variations in trigger time, a low photon count waveform has been modeled with different dark count rates. Each plot shows where a single trigger would occur for different random values draw in the Monte-Carlo model of the detection.

The first simulation shows where 5 possible triggers would occur (for random values of 0.10, 0.25, 0.50, 0.75 and 0.90) with a DCR of 0.

The second simulation employs the same random values, but with a DCR of 50 kHz. Note the bias in the cumulative probability in the lower time bins due to the increased noise. Even at this noise level, the location for these 5 possible triggers have not changed much compared to the previous no noise case.

The third simulation employs the same random values, but with a DCR of 250 kHz. Note the steeper slope in the cumulative probability in the lower time bins due to the increased noise. At this noise level, the random value of 0.10 has triggered outside (lower time bin) the plotted range and the triggers within the range have shifted to shorter time bins. Clearly the return near time bin #1150 would be considered a noise return since it does not correspond to any of the signal in the original waveform.

The GmAPD detector model is selected from the Detector Type menu on the Detector tab in the Lidar Detector Model. This provides the interface for the user to set the various parameters describing the detector’s characteristics.

The photo-multiplier tube (PMT) model emulates the detection system employed in the ATLAS sensor that will be on the ICESat-2 satellite. The ATLAS receiver uses commercial photo-multiplier tubes that feature a 4x4 segmented anode. The ATLAS system has 2 sets of PMT devices, with one set receiving the "strong beam" spots and the other set receiving the "weak beam" spots. The "strong beam" PMTs route each of the 4x4 anodes to a separate discriminator, resulting in 16 channels. The "weak beam" PMTs have groups of 2x2 anodes summed and then routed to separate discriminators, resulting in 4 channels. Note that the grid of anodes do not provide an imaging capability. A complete description of the ATLAS system can be found in external documents written by the system architect Anthony J. Martino at NASA’s Goddard Space Flight Center (GSFC).

Because of the photon counting nature of this PMT-based detection architecture, the implementation of this detector model is very similar to that employed in the Geiger-Mode detection scheme described earlier. The PMT detector system is characterized using the same parameters employed by the Geiger-mode scheme, including the photon detection efficiency (PDE) and the dark count rate (DCR). The major additions to the Geiger-mode model found in this model are the following:

The probability of a trigger in a given time bin i can be computed using the following equation:

\( P(i) = \left [ 1 -\exp(-S(i) -w) \right ] \)

where S is the received active signal (in photoelectrons) and w is the noise rate per bin (photoelectron events per bin from passive background flux and dark count). When comparing this equation to the Geiger-mode detection equation presented earlier, the Geiger-mode equation has an explicit trigger exclusion term (driven by noise events) that extends all the way back to the start of the listening window. For this model, the "dead time" (which is a trigger exclusion mechanism) is part of the Monte-Carlo simulation (described below). The support for multiple channels is implemented by simply dividing the received waveform N times (16 times for the "strong beam" PMTs and 4 times for the "weak beam" PMTs) and repeating the detection process on the weighted waveform N times.

The modeling technique is as follows:

The impact of dead time is fewer overall detections since the device will miss trigger opportunities during dead time. The decrease in the number of detections can be observed when comparing the "Dead Time = 0 ns" to the other dead times in the plot below. In addition, the peak trigger time statistics will start to skew toward "first photon" events because later events will be blocked by the dead time. The plot below shows this effect by plotting the normalized number of returns per time bin for a 2 ns wide return pulse.

The impact of multiple-channels is to increase the total number of detections and to compensate for the "first photon" bias introduced by "dead time" in strong return scenarios. The plot below illustrates how more channels not only increases the total number of triggers (the area under the curve) but also how additional channels help to mitigate the "first photon" bias (the peak shifts back toward 0).

The PMT detector model is selected from the Detector Type menu on the Detector tab in the Lidar Detector Model. This provides the interface for the user to set the various parameters describing the detector’s characteristics.

One of the more interesting areas of development right now for laser radar are "waveform" LIDAR systems. Unlike Linear-mode (Lm) or Geiger-mode (Gm) system that collect a relatively small number of height measurements per GSD, a "waveform" LIDAR (wLIDAR) essentially digitizes the returned flux during the time window at some time resolution. That means you have intensity information at all the ranges rather than just a discrete set of them. This allows for some very complex post-processing of the data. Most wLIDAR systems feature a high performance, single-pixel detector coupled to a high performance digitizer. To obtain area coverage that detection system is coupled to a scanner that can provide across-track scanning while the airborne platform provides the along-track scanning via platform motion.

The BIN file created in "Stage 1" by DIRSIG contains for each pulse, the temporal flux for each pixel in the receiver array during the period that the array is armed (defined by a range gate or listening window) and is accompanied by a lot of meta-data including the precise absolute time of the current pulse, ephemeris like the platform location, platform orientation, platform relative pointing and information about the source and receiver. That portion of the data containing the temporal flux is basically what a wLIDAR system collects. The major difference is that the data produced by DIRSIG does not include an impulse response of the detector, A/D system, etc. Unlike the point clouds generated by conventional systems that can use the standardized LAS format, there has not yet been a widespread adoption of a file format for waveform LIDAR data. Note that the LAS 1.3 specification released in 2009 indeed has support for waveform data, but most operating sensors are still using home-grown data formats that were created before the recent LAS additions. For example, the data could be formatted into an ENVI image file employing the "spectral" dimension as a temporal dimension. Additional files can carry important meta data for the pulse including location of the area images, the source pulse profile, etc.

Is there any easy way to convert a DIRSIG BIN file into an ENVI image file so that the waveform can be visualized and processed? The answer is "Yes" and it doesn’t require conversion. The BIN file is a 3D cube of data with some metadata. The ENVI image header file can be used to directly import the BIN as long as:

The uncompressed data requirement is because ENVI cannot read the compressed data. Pulse compression can be turned off in the Output tab of the Receiver configuration. The single pulse requirement is because ENVI cannot handle skipping over the meta data between pulses in single, multi-pulse file. If you have a multi-pulse simulation, you can tell DIRSIG to generate a separate file for each pulse using the Generate files per capture option on the Output tab of the Receiver configuration.

The required ENVI header file can be created by hand using the following template:

Here is an explanation of the various ENVI header file variables:

In order to go from Level0 time-of-flight data to Level-1 3D point clouds, each time-of-flight measurement must be converted to a range and then the platform location, platform orientation and platform-relative pointing must be used to locate that point in 3D space. The user can incorporate some common sources of "knowledge error" into this Level-0 to Level-1 processing.

The Platform and Mount tabs allow the user to incorporate errors into the perfect knowledge of the platform location, platform orientation and platform-relative pointing that is provided via the BIN file. The user can selective specify temporally uncorrelated standard deviations about the true values for the six degrees of freedom for the platform ephemeris and the three degrees of freedom for the platform-relative pointing.

The options for the ASCII/Text output include:

The ASCII/Text point cloud file writes a single line for each point in the point cloud. Each line (record) contains columns that describe the point. Depending on the detector model chosen and the options selected, the number of columns and what each column is changes. If the comment/header is included at the top of the file (the line begins with #), it describes the column (fields) for each point record (each line in the file). To clarify the format, a series of example files are provided below.

The LAS point cloud files contain the minimal headers and point meta-data descriptions described in the LAS 1.2 specification published by ASPRS.

An ASCII/Text point cloud file has the flexibility of being imported into most 3D plotting and visualization tools. The following is a partial list of the tools commonly available for basic visualization of ASCII/Text files:

An LAS format point cloud file can be loaded into large number of commercial LIDAR exploitation packages. Many of these software packages offer free versions with limited analysis capabilities. However, these free versions usually support basic loading and visualization of LAS format files. The following is a partial list of the free tools available for basic visualization of LAS files:

The LidarStatic1 demo shoots a single pulse onto a tree and receives that pulse onto a large format (128 x 128) array.

The LidarWhisk1 demo shoots a series of pulses while scanning in the across-track dimension.

The LidarDynamicGate1 demo explains how the range gate can be adjusted on a pulse-by-pulse basis.

The LidarUserPulse1 demo explains how a user-defined temporal pulse profile can be incorporated into a simulation.

The LidarPointing1 demo demonstrates how command errors (jitter) and knowledge errors (GPS/INS noise) impact 3D point clouds.

Under the hood, the Gaussian form utilized in the DIRSIG model is based on the form normally used for the Normal Distribution:

\( f(x) = \frac{1}{\sigma \sqrt{2 \pi}} e^{-\frac{1}{2} \left ( \frac{x}{\sigma} \right )^2} \)

where \(\sigma\) is the standard deviation.