LIDAR processing notes

Notes from visit to Leica, Switzerland on 11/Aug/2008.

Contact email:

  • sensors.support at the following host leica-geosystems.com
  • Various support engineers at Leica. The one that trained us was Barbara Haebler.

FTP site (for updates, etc):

Overall system design and comments of note

System is a Leica ALS50 (phase II) LIDAR. There is an accompanying 39 megapixel [7216x5412, 12bit] digital camera, referred to as the "RCD".

A GPS receiver and IMU (brand?) are included. The LIDAR, RCD and IMU, are mounted on a shock plate to protect it from strong movements - this will also isolate the IMU from the hyperspectral imaging sensors.

The IPAS controller (name correct?) allows for event recording, amongst other things. Phil plans to use this to record frame sync pulses from Eagle and Hawk.

There are also some associated control and display devices for operator and pilot usage.

How it works

The LIDAR works by firing a (4ns or 9ns) laser pulse downwards and measuring the roundtrip time for the light pulse to return, then converting this to a distance. The pulse isn't modulated by a carrier - it's just an on/off pulse. There are four timing cards ("range cards R1-R4") running for a pulse, so up to 4 returns can be detected (R4 actually detects the last return rather than the 4th?). The system has a "MPiA" (Multiple Pulses in the Air) mode, which fires two pulses evenly separated, rather than waiting for the first to come back before firing another [SPiA mode, times out in case the pulse is eaten]. To measure this, there are actually two banks of timing cards (bank A and bank B, both with R1-R4 cards), so there are 8 timing cards in total.

A minimum time separation between two returns means the minimum distance between two returns must be at least 2.7m for them to be counted as independent. The expectation for the number of returns is 1 return ~100%, 2 returns ~10%, 3 returns ~1%, 4 returns ~0.1% of points - obviously this varies with the terrain. When there are 4 returns, each range card measures the time of the return pulse. When there are less than 4 returns, R4 is a second measurement (not a copy of) of the last pulse - i.e. if there are 2 returns, you will have R1, R2 and R4 (= re-measurement of R2).

The intensity of a return is measured only for the first 3 returns (R1-R3), and is an 8 bit value (0=dark (water), 255=bright) relating to the reflectivity of the illuminated surface. The value is amplified by an automatic gain controller, and is not related to a physical measure (can it be?). The intensity can be used in various processing algorithms to help distinguish transitions between surfaces. The AGC tries to keep the intensity in the range 100-150 or so.

The laser is scanned across a (up-to) 75 degree swath by an oscillating mirror. Due to the acceleration/deceleration of the mirror, this produces a sinusoidal pattern to the trace on the ground, with the highest density of points at the peak and trough of the sine wave (i.e. at the edges of the swath). If the swath width is set to less than 75 degrees (45 degrees recommended), there's a roll compensation mechanism that tries to smooth out small roll movements by using the remaining freedom of motion. The scan rate of the mirror is dependent on the FOV angle (36Hz for 40 degree FOV, 24Hz for a 75 degree FOV, calculate with 412.33 * FOV-0.6548).

The laser is an 8W class 4 laser, operating in the infrared range. The divergence results in approximately a 22cm spot on the ground when fired from 1000m up.

Controllable parameters:

  • laser intensity (0 -> 8W output), controlled by operator as a percentage output. Has safety cutouts if the light level at the ground could cause eye damage.
  • altitude (kinda a parameter ;) ) - minmium of ~650 up to ~2000m (after 2km, you start getting poor returns on forests, etc, the real limit is up to about 6km in ideal conditions)
  • ... pulse frequency, scan angle, etc [TBD]

Automatic Gain Control

The intensity of a return is measured (R1-R3) as an 8 bit value (0=dark (water), 255=bright) relating to the reflectivity of the illuminated surface. The value is amplified by an automatic gain controller to keep it within 8 bit range. The value of the AGC is also measured as an 8 bit value and is recorded per pulse as the "gain control voltage" in the raw files. Typical values are 110-180.

The AGC is controlled automatically in our LIDAR, but operates on a threshold basis. If there are N (64?) low intensity points in a row, the AGC will step up one level. If there are M (32?) overbright points in a row, the AGC will step down. This may lead to visible steps in the image intensity, but isn't supposed to be a big problem. Look out for streaky areas though ("contact support if you see this" - though not sure why).

RCD camera

The system has a 39 megapixel digital camera, referred to as the "RCD", with the following characteristics:

  • 7216x5412 resolution, 12bit [=~57MB raw]
  • 60mm lenses (changeable), with 44.2 degree x 34 degree FOV
  • pixel size 6.8microns (= ~15cm ground resolution at 1350m)
  • 1/4000 exposure time, 2.2s per frame
  • serial number 21
  • shutter is curtain style, taking 8ms to open and close, so there's a 4ms time difference between the top and bottom of the CCD stopping acquiring light

Logfiles contain various handy info, notably the image event file contains GPS time, pos/att and picture id

Calibration files are required and supplied by Leica (dark, gain and camera parameters).

File name convention is: DDHHMMSSXXXXXXGN.raw

  • DD = day of month
  • HHMMSS = timestamp
  • XXXXXX = incrementing picture id
  • G = gain (1-4)

Orthophotos

This apparently takes another week long training course. We do have the basic software to do it but limited knowledge. Basic process is to determine the camera parameters, correct and overlay the images, keeping ideally only the most central parts. You then have to manually check through, especially along seams, and manually move the seam line in some instances (e.g. when perspective effects result in odd errors due to the seam hitting a tall building).

Recommendations

Perform calibrations frequently at the start of the period to get a feel for how they hold. Reduce down to fit circumstances over time.

Mission planning considerations

Water absorbs IR, so expect poor returns from wet surfaces. Ideally one should wait for a whole dry day after rain.

For differential GPS, aim to have a base station nearby (within 20km) and ideally in center of the scene. Base station data should be recorded at 2Hz

Ideally choose a time with a good GPS constellation (lots of satellites, PDOP >= 4).

Consider terrain and minimum target size to determine required point density. Steep terrain may cause shadowing effects due to perspective.

Consider reflectivity of surface - lighter surfaces need less laser power. Too high a power means intensity overflows on the reflection, which mess up the range - AGC should deal with this for large areas but it should be considered in planning, particularly for small bright targets on a generally darker background. Assume 10% reflectivity is typical.

Recommended to do a figure 8 loop at the start and end of the data acquisitions to make the IMU happy. Need one at the end for reverse navigation processing.

Processing

Required items for processing

Items required from every flight:

  • GPS and IMU data (including a basestation if not using PPP)
  • Raw laser data
  • Logfiles from LIDAR
    • Flight logfiles are useful too, if the LIDAR flight planning software is used
  • RCD & webcam images
  • Pressure and temperature measurements at the plane position above the site [this affects how long it takes the laser light to move through the air]

Items required in general:

  • calibration (see elsewhere)
  • lever arm measurements

Extraction

Use IPAS Pro to extract GPS, IMU and laser data from the raw files.

  • probably worth having the real time navigation info (doesn't add a lot of processing time)
  • during extraction, view the listing and verify the lever arms are correct (IMU one should never change).
  • look out for data gaps (listed in log)

Most Leica systems are mounted "laser backwards" (cables will be at the rear if this is so) - ensure that z=180 in the IPAS Pro aircraft tab.

See other pages on navigation processing as they cover this already.

General processing

Mark covered this better

QC

Things that may go wrong

In flight

Flying too low - if the LIDAR detects that the laser power may be too high for eye safety, it will cut out the laser automatically (if you're a seagull looking up, bad luck, it only accounts for ground height).

Flying outside the "range gate" (acceptable ranges of distances) may cause similar effects. Going too high will tend to make the edges of the swath drop out first (due to path length).

Automatic gain control weirdness - the measured intensity is returned via an AGC which may step up or down depending on the returns from the ground. Measured return intensity should only be used as a guideline rather than a real measurement.

Laser power too high - if the intensity overflows (reaches 255), the intensity based range correction will probably be wrong. For example, a freshly asphalted road with bright white reflective (overflowing) strips may appear with the strips appearing to float up to 20cm above the road surface.

Water absorbs IR, so expect poor returns from wet surfaces. Ideally one should wait for a whole dry day after rain.

In processing

Most Leica systems are mounted "laser backwards" (cables will be at the rear if this is so) - ensure that z=180 in the IPAS Pro aircraft tab.

Streaky intensity images - "bad, contact support for guidance". These will probably cause range errors due to the intensity variation, no idea what the cause would be as yet..

Poor accuracy at edges - consider trimming the swath width (reduce processed angle) to cut off bad edges.

Other error sources

Atmospheric effects; the light path may bend due to atmospheric refraction, density effects, etc causing the laser to hit the ground earlier or later than expected, or in a different position. This will be most noticeable at the edges of the swath where there are longer path lengths (and more atmosphere to pass through) and may look like the swath curls up or down at the edges (path length error), or may look like compression or stretching of the edges (if there is horizontal divergence). The error was claimed to be a fraction of a meter at 6km altitude (i.e. not well bounded). Measurement of temperature may help with this effect, but was said to be a minor value.

Range correction error; if the range correction is wrong, the electronics will measure the path length incorrectly. Points at nadir will be in error vertically only and points at the edge will have a vertical and horizontal error. The error will make a flat piece of land look like a smile and the nadir point will be below the expected land surface (one can only get "late" measurements?).

Torsion (of mirror) error; the mirror may be out of the expected position at the edges of the swath due to it bending under high acceleration. There will be no error at nadir (no acceleration) and errors at the edges, inducing a smile effect again, but with the nadir at the correct height. Range errors should be corrected before working on the torsion error, and the nadir point should be used for the range correction as there will be no mirror torsion effect there.

MPiA mode errors: if a seagull gets in the way of the second pulse before the first pulse has returned, things will mess up. On an edge of a very unluckily placed cloud, this would look a bit like the cloud merging into the ground. Presumably rare or minor.

Initial QC

See if the lines are too short or if the point cloud has poor return %s. If so, check the following (look at webcam images for hints):

  • clouds
  • height problem ("range gate" issue)
    • eyesafe shutoff (too close to ground)
    • too high may give dropouts
    • find altitude over ground (measured GPS alt - geoid-spheroid height) and see if it's within the min & max ranges

Laser power too high - if the intensity overflows (reaches 255), the intensity based range correction will probably be wrong. For example, a freshly asphalted road with bright white reflective (overflowing) strips may appear with the strips appearing to float up to 20cm above the road surface.

  • filter out incorrect points (based on high intensity and height difference from locality?)
  • also, if you observe odd spikes, check for intensity over 150 - this may indicate a two-peaked response with the first peak being 150 (and being the detected one) and the second being 255 (overflow), causing the effect above

Full QC procedures

Mark seems to have covered these better, and this will need some detailed clicky info, key concepts only here

Essentially all using TerraScan. The main things are:

  • profiles
    • profiles along averaged areas
    • profiles along a track (e.g. following a ditch / feature in the surface)
    • narrow profiles across features where swaths overlap
  • comparisons against GCPs
    • requires the point cloud to be converted to a surface (

Calibration

Parameters needing calibration

Boresight parameters (pitch, roll, heading).

  • Angle between straight-down and what the sensor thinks is straight-down, as it's mounted in the plane.
    • Angle between a line from the sensor head (mirror centre) to the point on the ground at the centre of the swath and a line from the sensor head to the centre of the spheroid (or the reference frame's Z axis?).
  • Measured in calibration procedure
    • We need to do this

Pitch error slope.

  • The mirror will not be mounted exactly flat to the laser so, as the mirror moves, the pitch of the beam will change by a small amount.
    • Roll and yaw either have no error slope (laser position central? geometry means no effect?) or a negligible effect (presumably, as there's no parameter)
  • Measure pitch at nadir and at swath edges to determine how the pitch changes - the first order for this is is the pitch error slope.
  • Measured in calibration procedure
    • Correct pitch error (at nadir) first
    • We need to do this

Lever arm measurements:

  • IMU -> sensor head (centre of mirror?)
    • Measured in the factory and provided by Leica
  • IMU -> measurement point [front left corner of the casing]
    • Measured in factory and provided by Leica
  • measurement point -> GPS antenna centre
    • We need to measure this one with a survey

Range offset correction (+range card calibration).

  • Correction for the slightly different timing of the 4 range cards in both banks in the system.
  • At a set distance, the range cards should all return the same result.
  • Measured by Leica but also measured and verified in calibration procedure. (see below)

Mirror torsion parameter.

  • This measures how much the mirror flexes at the edges of the swath, when under high acceleration, putting it out of position.
  • Later Leica mirrors (including ours) are much stiffer and apparently do not flex.
  • Fixed very high value (100,000 N?) meaning "no measurable flex".
  • Can be measured during the calibration procedure, but should never be necessary.

Factory calibrated

Intensity Based Range Correction (IBRC).

  • A range correction based on the reflectivity of the target object ("reflections are slower on darker objects"? perhaps this is more that a lower intensity return may take longer to be detected - i.e. takes longer for enough photons to come back to exceed the detection threshold).
  • Measured in factory and provided by Leica.

Transitional Pulse Rate.

  • 2 pulse widths are used (4ns and 9ns) depending on the pulse frequency (higher frequencies require a shorter pulse time), with the system switching over to a shorter pulse at a set frequency (around 120KHz?). This parameter measures the time difference for the pulses' return.
  • In our system, this distance is of the order of 11cm.
  • Measured in factory and provided by Leica.

Encoder Offset (scan angle correction).

  • The rotating mirror has a (optical?) shaft encoder attached that tells it the position of the mirror. The mirror will not be mounted exactly at nadir relative to the shaft encode and this value measures the offset between the centre of the scan pattern and the centre point of the encoder.
  • The value is measured in "ticks" and will appear quite high as the shaft encoder has a high resolution (order of millions of ticks).
  • Measured in factory and provided by Leica.

Calibration site requirements

Flight pattern is a cross made of 4 opposing flight lines in the 2 directions (e.g. N->S v S->N + W->E v E->W), and one set of parallel lines with 50%(?) overlap.

The features required for a perfect cal site are:

  • Need a source of multiple returns - tallish (15m) trees in a forest are best. Try to include a treed/forested area in some parts of the flight lines (doesn't have to be in all, nor in the central area). The multiple returns are required for the range card calibration (see below).
  • Straight, flat areas made of a hard substance that'll generate only one return pulse. Roads or runways ;) These are used in the boresight (roll) and range calibrations.
  • Sloping areas are used to detect pitch and heading errors in the boresight calibration.
    • A road running up a hill is good for pitch (vertical change in surface easily found in the image). These are good for measuring along-track shifts.
    • Slopes with peaks (house rooftops are ideal) with the peak cutting across the line of flight ( ---> /\ ). Again good for measuring along-track shifts.
  • An accurate ground survey (see below).
  • Slow overflight for maximum point density.

For range offset calculation

GCPs required, ideally 1cm vertical accuracy (between all measurements rather than absolute?).

Around 30-40 GCPs should be within the area covered by a 14 degree swath (+/- 7 degrees either side of nadir) taken at an altitude of ~750m [= ground swath width of about 180m].

General accuracy measurements

Other GCPs (number?) should be scattered around a wider area within the full swath width - typically run at 45 degrees [=620m wide on ground @ 750m alt] or to a max of 75 degrees [=1150m], though 75 degrees will introduce more errors..

Quality, accuracy, etc

A nominal quality of 5-10cm (vertical) is suggested as a reasonable output. If one has GCPs, it should be possible to do better (EA claims 1-2cm).

Range offset calibration

Range offset correction (+range card calibration) is to correct for the slightly different timing of the 4 range cards (R1-R4) in both banks (A & B) in the system, and to correct any overall ground offset.

Two datasets are required:

  1. A real dataset with:
    1. a source of multiple returns, such a forest (for step 1 below)
    2. a strip with well known distances (for verification)
  2. BIT (Built-In Test) mode data, where the range cards are all electronically fed with identical fake data representing the same distance. All cards should give the same result, so differences are used to calibrate each card against the others.

First, we need to determine the timing differences between the 4 range cards (R1-R4) in each bank. To do this, we use a dataset with multiple varying returns present - we need combinations of 2-4 returns (forests are good for this, being tall and porous enough to give multiple returns). When there are less than 4 returns, R4 is a second measurement (not a copy of) of the last pulse - i.e. if there are 2 returns, you will have R1, R2 and R4 (= re-measurement of R2). This is exploited to compute the difference in timing between R2 and R4 (averaging many 2-returns). Similary 1-return and 3-return pulses are used to measure R1-R4 and R3-R4 differences. The end result is a set of timing differences between all the cards in a bank.

Second, we need to establish the timing differences between bank A and bank B. As with the first step, we need the timing cards to measure exactly the same instant. To do this, we use BIT mode data, where R1 in both bank receives the same electronically generated pulse at the same instant and are thus measuring the same event. Averaging these numbers gives the timing offsets between the R1 cards, which can be combined with the first measurements to establish timing between all cards.

Procedure:

  • run the timing estimate program on the two datasets above, using a nominal range error of 0
    • hopefully Mark wrote down the details on how to do this - it's somewhere off the ALS preprocessor menus, but looked pretty straightforward. Resulting numbers need to transcribed to ALS processor afterwards
    • do we need to turn off "use average of last returns?"
  • process a 14 degree strip (+/-7 degrees of nadir) over the dense GCP region (need ~30-40 GCPs @ 1cm vertical accuracy)
    • in TerraScan, compute the average error between the GCPs and surface generated from the point cloud (didn't write this procedure down either)
  • enter this error as the new nominal range error into the timing estimate program and re-run
  • reprocess the strip and verify against the GCPs - error should be dz = ~0 (<1cm), stddev ~5cm.

Boresight calibration

The basic method is similar to the imaging instruments, but exploits the measurement of range. A set of flight lines and data points are chosen, in a particular order, such that specific calibration parameters are isolated.

The flight pattern for boresight calibration is a cross made of 4 opposing flight lines in the 2 directions (e.g. N->S v S->N + W->E v E->W), and one set of parallel lines with 50%(?) overlap.

Manual procedure and QC of automatic procedure

Pitch

Take opposing flight lines. If there is a pitch error, the lidar will be sensing pulses ahead/behind (along track) the nadir point. There will be a small overall range error, but the main effect will be that the data is shifted backwards/forwards of where it should be.

Search for areas close to nadir where the height of the surface changes in an clear and predictable fashion, with reasonable sampling density (e.g. a road leading up a hill, or rooftops with their peaks going across-track).

Take along-track profiles on both opposing flightlines.

  • If there is no pitch error, there will be no offset and the resulting profiles will coincide (the hill profiles will overlay and look like /).
  • If there is a pitch error, the profiles will be offset along track and the resulting profiles will not coincide (the hill profiles won't overlay and will look like / /). By measuring the along-track offset, one can estimate the pitch error.

One must use profiles close to nadir to avoid being affected by the pitch-error-slope parameter (see below).

Repeat until the pitch error has been removed.

Roll

If there is a small pitch error, surfaces will be tilted up on one side of the swath and down on the other. In opposing lines, these errors will be on different sides.

  • Take opposing flight lines.
  • Find flat surfaces going across the swath (e.g. a road going across-track).
  • Take across-track profiles on the flat surfaces
  • If there is a roll error, the profiles will not align with one another and will have a X shape rather than an overlapping - shape (the flatter the X, the smaller the error)

Vertical displacements due to roll error will be most noticeable at the edges of the swath. For small errors, the nadir point will be roughly the same, for large ones, there will be some across-track shift, but that will be much harder to measure than the vertical errors.

Heading / Yaw

If there is a heading error, the left side of swath will move backwards while the right side will move forwards (or vice versa). If you reverse the direction of flight (ie. with opposing flight lines), these errors will overlap (left side down in one direction matches right side up in the other) and not be visible.

  • Instead, take two parallel flightlines with 30-50% overlap (more overlap = less chance of detection)
  • Along the mid line of the intersection, find areas where the height of the surface changes in an clear and predictable fashion, with reasonable sampling density (e.g. a road leading up a hill, or rooftops with their peaks going across-track).
  • If there is an along-track shift between features, there is a roll error (size of shift proportional to error)
  • If there is no error, features should overlap

Pitch error slope

The mirror will not be mounted exactly flat to the laser so, as the mirror moves, the pitch of the beam will change by a small amount. The resulting errors are along-track offsets just like the pitch errors, but varying, with no error at nadir and maximum error at the edges of the swath.

To correct this, first correct the pitch error at nadir (the overall pitch error), then determine the pitch error at the edges of the swath. This will let you compute an error slope (first order linear).

  • Straight, flat areas made of a hard substance that'll generate only one return pulse. Roads or runways ;) These are used in the boresight (roll) and range calibrations.
  • Sloping areas are used to detect pitch and heading errors in the boresight calibration.
    • A road running up a hill is good for pitch (vertical change in surface easily found in the image). These are good for measuring along-track shifts.
    • Slopes with peaks (house rooftops are ideal) with the peak cutting across the line of flight ( ---> /\ ). Again good for measuring along-track shifts.

Mirror torsion error correction

This shouldn't need to be done with newer LIDARs, but included here in case. Leica suggests checking for torsion errors if there seem to be problems in the boresight calibration and reporting to them if found.

Visible sign of torsion error is a smile effect.

Detect torsion by:

  • choose highest altitude pair of crossing lines
  • pick a point at nadir on one swath and at the edge on the other
    • torsion will be at max for the edge swath and at 0 for the nadir swath
  • take cross sections to see if there is any consistent vertical difference

Automatic method (using Attune)

Needs per click instructions as it's a bit clunky - I have some of these, but hopefully Mark has more

Note that Attune apparently crashes if you use >400MB files, and crashes plenty anyway - save often! Do an initial run with just a few tie points to make sure the lines will work.

First, classify the LAS files - the automatic process should only be used on ground points (largely to avoid perspective effects on buildings?)

Load the LAS files into Attune and georectify, picking an appropriate pixel size.

  1. Select ~50 - 100 tie points (see below)
  2. Enter any pre-known/estimate boresight parameters (mirror torsion 100000N) to give it a head start.
  3. Enter image observation weights (x=0.2, y=0.2, z=0.1m), indicating how accurate we are in each dimension
  4. Select "use class 2" (ground points), max iterations=50, angular convergence 1x10E-08
  5. "Solve" (don't compute torsion parameter)
  6. Examine aposteriori reference -> if it's ~1-2, stop, otherwise you need to fix bad tie points
    1. Aim for a standard deviation of ~0.00001 rad (P & R) and ~0.00005 (H)
    2. To fix bad tie points, sort by the residual error and delete or correct the worst points
    3. A bad tie point will have a residual error of > 25cm (X & Y) or > 5cm (Z)
    4. Rerun and try again..

After completing the Attune process, save the parameters as a reg file (for reference - you will need to manually transfer the numbers to avoid overwriting other parameters). Then check the boresight parameters using the manual procedure above and adjust if required.

Selecting tie points

  • Better to have fewer good tie points than many bad ones
  • Focus mainly on the lower altitude lines
    • Place tie points on the higher altitude lines only when you're really sure of them - they will always have a greater error due to the pixel size.
  • Tarmac roads are good, but:
    • avoid paint on road (intensity problems and height difference)
    • choose continuous homogeneous areas
  • Avoid:
    • building corners (perspective)
    • grass (multiple returns, moves)
    • cars/people (moves)
    • areas with streaky intensity (likely to have ranging errors)
  • Try to select areas where you can use the nearby geometry to help you
    • e.g. follow straight road edges to where they intersect another feature, so you have more than just one pixel to guide you to the intersection point in other images
    • avoid using building edges (perspective shifts)
  • Location of tie points
    • Focus mainly in the core intersection area
    • Points around nadir
    • Points around the edges of the swath
    • Some generally scattered around
    • Some outside the core intersection area
    • Try to include some on slopes(?)

Ground classification

Classification is the process of labelling the points in the point cloud by type. Ground classification specifically refers to labelling points that're considered to be part of the ground(!) This is a required step for the automatic boresight procedure.

Steps:

  • remove very low, noisy points from consideration first
    • e.g. remove groups of up to 3 points > 1m below neighbours (iterate a few times so it settles)
  • classify groups of low points as ground (run N times to settle), takes account of neighbours
  • repeat with singletons
  • remove isolated points (>5m from any other point)
  • manually "add points to ground" for areas that should be ground but haven't been classified as such (e.g. peaks of hills) and rerun process above

Run as a macro. Terra->tools->macro->new.

Random snippets

SCN file, DC=delta counter (time in ms since last GPS-second tick), ANG=angle in ticks, RI = return? (in m?).

Last modified 15 years ago Last modified on Aug 28, 2009, 10:47:20 AM