Automatic Ice Thickness Estimation from Polar Subsurface Radar Imagery

1
Automatic Ice Thickness Estimation from Polar
Subsurface Radar Imagery

• Gladys Finyom
• Michael Jefferson Jr.
• MyAsia Reid
• Christopher M. Gifford
• Eric L. Akers
• Arvin Agah

2
Overview

• Introduction
• Background/ Related Works
• Overview of Remote Sensing
• Challenges of Processing Radar Imagery
• Ice Thickness Estimation from Radar Data
• Methods
• Edge Detection and Following Approach
• Active Contour Cost Minimization Approach
• Experimental Results
• Conclusion

3
Introduction

• Remote sensing methods
• CReSIS uses Radar and Seismic/acoustic to acquire
  subsurface data from a remote location. (i.e.,
  surface, air, or space).
• Radar and Acoustic sensors
• used to gather data about the internal and bottom
  layers of ice sheets, from the surface.
• Other Examples
• Satellite-based imagery, and identification of
  events.

4
Introduction

• Surfaced-based and Airborne radio echo sounding
  of Greenland and Antarctica ice sheets
• Determine ice sheets thickness
• Bedrock Topography (smooth, rough)
• Mass Balance of large bodies of ice.
• Challenges in Radar Sounding
• Rough surface interface
• Stages of melting (top, inside)
• Variations of ice thickness, topography
• Processing Data
• Requires knowledge about sensing medium
• Ultimately used for scientific community

Greenland's ice sheet
NASA/Rob Simmon
5
Introduction

• Goal
• Focus on automating task of estimating ice
  thickness.
• Process
• Identifying and accurately selecting of ice
  sheets surface, interface between the ice, and
  the bedrock.
• Knowing the surface and bedrock in the radar
  images
• helps compute the ice thickness.
• help studies relating to the ice sheets, their
  volume, and how they contribute to climate change.

Four outlet Glaciers studied by CReSIS
researchers. Leigh Stearns
6
Overview of Radar Remote Sensing

• Radars transmit energy in form of a pulse from
  an antenna, energy reflects off of target(s), and
  is received by an antenna.
• Distance measured based on energy travel time
  back and forth from the targets.
• Gives reflection intensity and depth information
  about the targets.
• Ground Penetrating Radar (GPR) able to observe
  properties of subsurface, ranging from soil,
  rock, sand and ice.
• When data is collected, the targets are internal
  layering in the ice sheets which have a strong
  echo return from the bedrock beneath the ice.
• Interface (3.5 km or gt) below the surface.
• Requires great transmit power and sensitive
  receive equipment because of energy loss within
  ice and with depth.

7
Overview of Radar Remote Sensing

• Each measurement is called a radar trace, and
  consist of signals, representing energy due to
  time. The larger time correlates with deeper
  reflections.
• In an image, a trace is an entire column of
  pixels, each pixel represents a depth.
• Each row corresponds to a depth and time for a
  measurement, as the depth increases further down.

8
Overview of Radar Remote Sensing

• A flight segment consist of a collection of
  traces which represent all the columns of the
  image, from the beginning (left) to the end
  (right) during flight.
• A pixel width represents the track distance
  between traces, and depend on the speed of the
  aircraft during the survey.
• The flight segment is called an Echogram.

9
Overview of Radar Remote Sensing

• The Energy from the radar into the ice changes in
  dielectric properties (air to ice, ice to bed
  rock) and causes the energy to reflect back.
• Water surrounded by the ice, and frozen ice
  against the bedrock both represent a strong
  reflecting interface.
• To determine whether each is present, it depends
  on the radar and its setting.

Reflection intensities are strongest at the
surface and weaker because of depth. Depth
increases from left to right.
10
Example Radar Echogram Greenland 05/28/2006
Figure shows radar echogram over an ice sheet,
illustrating the reflection of internal layers
and the bedrock interface beneath the ice sheet.
11
Challenges of Processing Radar Imagery

• Automated processing and extraction of high level
  information from radar imagery is challenging.
• Noise is usually electromagnetic interference
  from other onboard electronics.
• Low magnitude, faint, or non-existent bedrock
  reflections occur
• Specific radar settings
• Rough surface/bed topography
• Presence of water on top/internal to the ice
  sheet.
• This produces gaps in the bedrock reflection
  layer which must be connected to construct an
  adjacent layer for the completion of ice
  thickness estimation.

12
Challenges of Processing Radar Imagery

• Backscatter introduce clutter and contributes to
  regions of images incomprehensible.
• In addition, bed topography varies from trace to
  trace due to rough bedrock interfaces from
  extended flight segments.
• Lastly, a strong surface reflection can be
  repeated in an image, surface multiple due to
  energy reflecting off of the ice sheet surface
  and back again.
• If there is a time difference between the first
  and second surface return, the surface layer will
  repeat in an image, at a lower magnitude with an
  identical shape.

13
Ice Thickness Estimation from Radar

• Along with raw data values, and GPS location
  measurements, ice thickness is needed for
  scientist to
• study mass balance
• sea level rise
• Environmental/ human impacts
• Ice thickness is computed by selecting the
  surface and bedrocks reflections in pixel/depth
  coordinates, for each trace, and subtracting
  their corresponding depths.
• Several experts had select the layers using
  custom software.

14
Ice Thickness Estimation from Radar

• The surface is selected based on the first and
  largest reflection return.
• The bedrock, is more challenging due to buried
  noises.
• Experts tend to skip traces to speed up the
  process.
• This causes errors, and require more time to
  estimate ice thickness in a single file.
• Thousands of images need work, and the manual
  approach is not sufficient.

Figure shows CReSIS picking software, the surface
return is fully picked, while bedrock return is
partially picked.
15
Related Work

• Internal Layer
• predicts depth in certain layers.
• depth and thickness of Eemain Layer in Greenland
  ice sheet utilized Monte Carlo Inversion flow
  model to estimate unknown parameters guarded by
  internal layers. 1,2
• Identification
• Layers, contours, and curves are done using image
  processing, and computer methods.
• Adaptive contour snake fitting, where an image is
  a cost grid, and represents a certain amount of
  energy. 3,4
• Such approaches have been used in medical imagery
  (such as MRIs and CAT scans). 5

3 T. F. Chan, L. A. Vese, Active Contours
Without Edges, IEEE Transactions on Image
Processing 10 (2) (2001) 266277. 4 M. Kass, A.
Witkin, D. Terzopoulos, Snakes Active Contour
Models, International Journal of Computer Vision
(1988) 321331. 5 J. Kratky, J. Kybic,
Three-Dimensional Segmentation of Bones from
CT and MRI using Fast Level Sets, in Medical
Imaging Proceedings of the SPIE, Vol. 6914,
2008, pp. 69144769144710.
1 S. L. Buchardt, D. Dahl-Jensen, Predicting
the Depth and Thickness of the Eemian Layer at
NEEM Using a Flow Model and Internal Layers,
in Geofysikdag, 2007. 2 S. L. Buchardt, D.
Dahl-Jensen, Estimating the Basal Melt Rate
at NorthGRIP using a Monte Carlo Technique,
Annals of Glaciology 45 (1) (2007) 137142.
16
Edge Detection and Following Approach

• Introduction
• Edge Detection, Thresholding, Edge Following
• Surface should be max value in each trace
• Bedrock should be the deepest contiguous layer in
  image
• Similar Work
• Skyline Detector
• Growing seeds in the sky
• Identify week clouds in sky imagery
• Our Approach
• Trace processed from bottom-up fashion until a
  strong edge is encountered

17
Automatic Surface and Bottom Layer Selection

• Surface Selection
• Extracting the location of the ice sheet surface
• The depth corresponding to the max value of each
  trace is selected as location of surface
  reflection
• Bottom Selection
• Preprocessed by
• Detrending
• Low-pass filtering
• Contrast adjustment

Figure Echogram that has been preprocessed using
detrending, low-pass filter, and contrast
enhancement
Figure Normalized echogram gradient magnitude,
showing the image edges
18
2D Derivative of Gaussian Kernel
Figure 2D derivative of Gaussian convolution
kernels (1.5 ) for computing vertical (left) and
horizontal (right) image gradients.
19
Cleaned Edge Image Result Image
Figure Echogram with overlaid automatically
selected surface (top, red) and bedrock (middle,
blue) layers using the edge-based method.
Figure Cleaned edge image following
thresholding, morphological closing and thinning
operations
20
Active Contour Cost Minimization Approach

• Similar work
• - Mars Exploration Rovers (MER)s automatic sky
  segmentation system
• - Further analysis of segmentation
• Our Approach
• Contour technique to fit a contour to the bottom
  layer

21
Automatic Surface and Bottom Layer Selection

• Surface Selection
• Same as Edge Detection technique
• Bottom Selection
• - Data preprocessing
• - EdgeCosts 1/v(1Gradient Magnitude)?
• - Creating an Image Gradient for upward force
• - Adding the edge cost image and upward force
  image

Figure Edge cost image, enforcing low cost for
strong edges and high cost for noise regions
22

• Initialization
• The contour is allowed to adapt until it reaches
  equilibrium
• 2N1 window (N 50 pixels) is maintained
• Window utilized for computing local stiffness to
  instill continuity and smoothness during
  adjustment
• Determination of Lowest cost (0) pixels and
  Highest cost pixels (1)?
• This technique allows the contour to fit to the
  image and bridge gaps in the bedrock layer

Bottom Layer Selection Continues
Figure Combined edge cost and upward cost images
23
Contour Cost Window Active Contour
Figure contour stiffness cost window during
processing (left) for the contours configuration
during the 75th iteration (right), illustrating
how the contour is encouraged to make smooth
transitions from trace-to-trace.
24
Final Adjustments

• TotalCostWindow(t) EdgeCosts(s) a x
  UpwardCosts(t) ß x ContourStiffnessCosts(t)?
• The minimal pixel location at each trace is
  selected as the contours starting configuration
• If the configuration does not change between
  iterations, or 500 iterations have been
  processed, the contour is determined to have
  reached equilibrium
• Ice thickness is computed for each trace by
  converting pixels for the bedrock selection to a
  depth in meters and subtracting it from the
  surface depth for each trace

Figure Echogram with overlaid automatically
selected surface(top, red) and bedrock (middle,
blue) layers using the active contour method.
Green is the initial contour configuration.
25
Active Contour Configuration
Figure Example contour adaptation sequence
throughout processing, illustrating how the
contour adapts to the bedrock interface and fits
itself to the most salient edge near the bottom
of the image
26
Experimental Results

• Processed in Matlab
• Data were 15 random subsets of 75 extended
  flights from Greenland.
• Range from 800-3000 rows and 1750-14500 columns
  (traces).
• Previous manual selection method took roughly 45
  minutes per file with approx 7500 columns per
  file.
• Automated edge-based method takes 15 seconds per
  file.
• Active contour (snake) method takes 2.5 minutes
  per file.

27
Results

• We assumed that the human approach was 100
  accurate
• Selection is considered correct if it is within
  5 of the human selection
• There are drawbacks with the manual approach

28
Edge-Based Method

• This method differs slightly from the active
  contour results even though both used same
  technique
• No Continuity
• Active Contour Method
• Method is able to outperform the edge-based
  method
• Takes a little longer to process images
• Smooth
• Continuous

29
Edge Method vs. Contour Method
Gap in bedrock
Contour method bridges the gap
30
Continued Example
Plotted points above the bedrock
Active Contour method rids the echogram on
non-continuous plotted pixels
Plotted pixels below actual bedrock
31
Continued Example
Edge-detection method works better
Artifact/Noise in the bedrock layer
32
Human Expert vs. Edge Method vs. Active Contour
Method
33
QUESTIONS / COMMENTS
34
34