Wasklewicz et al. Progress in Earth and Planetary Science (2017) 4:32
DOI 10.1186/s40645-017-0144-7

RESEARCH ARTICLE

Progress in Earth andPlanetary Science

Open Access

Simulating and quantifying legacy
topographic data uncertainty: an initial stepto advancing topographic change analyses
Thad Wasklewicz1* , Zhen Zhu2 and Paul Gares1

Abstract
Rapid technological advances, sustained funding, and a greater recognition of the value of topographic data havehelped develop an increasing archive of topographic data sources. Advances in basic and applied research related
to Earth surface changes require researchers to integrate recent high-resolution topography (HRT) data with thelegacy datasets. Several technical challenges and data uncertainty issues persist to date when integrating legacy
datasets with more recent HRT data. The disparate data sources required to extend the topographic record back intime are often stored in formats that are not readily compatible with more recent HRT data. Legacy data may also
contain unknown error or unreported error that make accounting for data uncertainty difficult. There are also casesof known deficiencies in legacy datasets, which can significantly bias results. Finally, scientists are faced with thedaunting challenge of definitively deriving the extent to which a landform or landscape has or will continue tochange in response natural and/or anthropogenic processes. Here, we examine the question: how do we evaluate
and portray data uncertainty from the varied topographic legacy sources and combine this uncertainty with current
spatial data collection techniques to detect meaningful topographic changes? We view topographic uncertainty as
a stochastic process that takes into consideration spatial and temporal variations from a numerical simulation andphysical modeling experiment. The numerical simulation incorporates numerous topographic data sources typicallyfound across a range of legacy data to present high-resolution data, while the physical model focuses on more
recent HRT data acquisition techniques. Elevation uncertainties observed from anchor points in the digital terrain
models are modeled using “states” in a stochastic estimator. Stochastic estimators trace the temporal evolution of
the uncertainties and are natively capable of incorporating sensor measurements observed at various times inhistory. The geometric relationship between the anchor point and the sensor measurement can be approximatedvia spatial correlation even when a sensor does not directly observe an anchor point. Findings from a numericalsimulation indicate the estimated error coincides with the actual error using certain sensors (Kinematic GNSS, ALS,
TLS, and SfM-MVS). Data from 2D imagery and static GNSS did not perform as well at the time the sensor is integrated
into estimator largely as a result of the low density of data added from these sources. The estimator provides a historyof DEM estimation as well as the uncertainties and cross correlations observed on anchor points. Our work provides
preliminary evidence that our approach is valid for integrating legacy data with HRT and warrants further explorationand field validation.
Keywords: Stochastic models, Digital terrain models, Legacy data, Geomorphic change detection, High-resolution
topography* Correspondence: wasklewiczt@ecu.edu
1
Department of Geography, Planning, and Environment, East Carolina
University, Brewster Building A227, Greenville, NC 27858, USA
Full list of author information is available at the end of the article© The Author(s). 2017 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0
International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, andreproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to
the Creative Commons license, and indicate if changes were made.

Wasklewicz et al. Progress in Earth and Planetary Science (2017) 4:32

Introduction
Topographic data sources are an essential source of information used to not only address scientific questions
related to landscape evolution, but also to inform policy
makers and land managers of recent environmental
changes. Earth scientists have also begun to carve out arole in the global climate agenda (Lane 2013) and topographic data will play a significant role in these types ofanalyses. Recent examples of this research avenue place
greater emphasis on understanding how humans modify
the physical environment (Tarolli and Sofia 2016) as well
as examining the likelihood of extreme weather leading
to magnitude and frequency changes in geomorphic processes and surficial features (Spencer et al. 2017). These
two recent research trends also highlight a need to beable to examine recent topographic changes in the context of longer-term rates of change. This task will require greater reliance on integrating legacy and more
recent high-resolution topography (HRT) data sources
to aid in our understanding of environmental change(Glennie et al. 2014; Pelletier et al., 2015). However, theintegration of disparate topographic data sources introduces numerous opportunities to increase uncertainty inmeasurements through space and time. Herein, we provide a means to assay and portray data uncertainty whenfusing legacy and recent HRT data to assess change in
landforms. Statistic and stochastic models are employedthat are not only consistent with topographic data structures simulated in this research, but also can track thecorrelation of errors over time and across different data
sources to establish reliable detection levels. Our findings represent an initial step to providing researcherswith the ability to detect topographic changes with a degree of certainty.
Topographic legacy datasets come from a variety offield and remotely sensed sampling techniques and instrumentation (Wasklewicz et al. 2013). Data quality,
resolution, and temporal availability from these disparate
sources have varied significantly over their historical development. Current research challenges in earth science
rely on a deeper understanding of how to integrate thesevaried datasets and their associated data uncertainty.
Advancement of scientific research on landscape evolution and environmental change detection rely on definitively measuring how the changing earth surface. Moreexacting measures would aid in items such as the development of more precise early warning systems wherelandscape dynamics pose hazards and risks to society,and could lead to potential innovations in the design ofinfrastructure that is more resilient to dynamic surficial
processes.
Research focused on topographic changes has reliedheavily on pre-event topographic data sets (e.g.,
Wasklewicz and Hattanji 2009; Wheaton et al. 2010).

Page 2 of 13

Some studies use the same instrumentation and apply
a consistent methodology throughout the observation
period (e.g., Wester et al. 2014; Staley et al. 2014;
Wasklewicz and Scheinert 2015). The integration of
topographic data from a single source with the aid ofa repeatable surveying campaign over time has presented an opportunity to reduce the systematic errors
while accounting for other positional errors and surface representation uncertainties. However, when repeatable surveying campaigns are not followed, researchers
found data can possess erroneous calibration and improper error modeling (Oskin et al. 2012; Glennie et al.
2014). These inherent issues are expressed as substantial
systematic errors, which lead to improper measurements
when compared with the post-event data. Glennie et al.
(2014) warn systematic errors must be minimized or removed prior to differencing the pre- and post-event airborne laser scanning (ALS) data sets. Schaffrath et al.
(2015) identified comparable issues with both verticaland horizontal measures from pre- and post-flood ALS
data being inadequately defined from the use of different geoid models and poor co-registration of flightlines, respectively.
The complexity of the data uncertainty increases with
the incorporation of legacy data. A single instrument isnever used in these cases. Instead, researchers attempt
to fuse topographic data from multiple instruments astechnological advances in data collection have occurred
over time (Crowell et al. 1991; Carley et al. 2012; Schaffrath et al. 2015). A variety of topographic data sourcesthat include contour maps, cross-sections or topographic profiles, raster, triangulated irregular networks(TINs), and point clouds can be integrated as legacy data
into these types of analyses. Each different data source
introduces varying: data quality, spatial resolution, andtemporal consistency. These items increase the spatiallyvariable uncertainty of the measurements taken during
analyses of these disparate data sets, which must be
accounted for in the presentation of the findings.
The recent incorporation of SfM-MVS (structure from
Motion-Multi-View Stereo) techniques has the benefit ofpermitting researchers to extend back further the legacy
record by using archived aerial photographs to reconstruct a point cloud and digital elevation model (DEM)(Gomez et al. 2015; Micheletti et al. 2015), but also adding further complications to examining spatially variable
uncertainty. While all signs from these studies indicate astrong potential to use these data sources to measure
topographic changes from legacy data, several issues
must be addressed to accurately use this approach. Someof the common photogrammetric issues included inconsistent image quality, varied scales, objects in motion,
clouds, and other superfluous information in the photos(Gomez et al. 2015; Micheletti et al. 2015). Initial results

Wasklewicz et al. Progress in Earth and Planetary Science (2017) 4:32

indicated aerial imagery at a scale of 1:20,000 may only
be able to detect changes in the 1 to 1.5 m range withquality ground control points (Micheletti et al. 2015),
but in some instances this value can be at a coarserscale (Fonstad et al. 2013; Gomez et al. 2015). Anothermajor concern with applying this technique is overestimation of the topography within areas of vegetation andlocations where there is disparate topographic relief(Gomez et al. 2015; Micheletti et al. 2015). However,some these obstacles can be overcome with the appropriate establishment of ground control points and consistent field validation of the results where possible andothers will require development of uncertainty measurescapable of uncovering errors in complex topography.
Robust estimates of spatially variable uncertainty havereceived more recent attention, but remain in their infancy (Schaffrath et al. 2015). As highlighted in Carley etal. (2012), researchers have generally used three different
approaches when examining uncertainty in topographic
legacy data: (1) a uniform error threshold (Brasington etal. 2003; Lane et al. 2003; Milan et al., 2007), which havebeen noted to bias data in various topographic sequences; (2) spatial variance thresholds used to produce
a minimum level of detection raster file (Heritage et al.
2009; Milan et al. 2011), which work well in settings
with high-density point clouds and mid- to fine-scaleroughness features; and (3) error-source thresholdmethods developed from fuzzy inference systems(Wheaton et al. 2010, Schaffrath et al. 2015), which has
been proven to assess the spatial variability of uncertainty often inherent in the digital elevation models(DEM[s]) used in topographic changes. Carley et al.
(2012) employed a hybrid method of the spatial varianceapproach (#2) to consider the addition of legacy topographic map data to their analyses. They applied a survey and interpolation error (SIE) equation (Heritage etal. 2009) to produce SIE point grids for the DEMs usedin the DEM-of-difference and combined these to produce
a critical error threshold (Brasington et al. 2000, 2003;
Lane et al. 2003) for each cell and a level-of-detection surface was generated from this style of analysis.
The issue of incorporating additional measurementswith existing models or maps has also been encounteredin the robotics and computer vision communities. For
example, a robotic mapping problem can be solved byusing discretized map representation. A common approach of mapping a 2D or 3D space is to abstract thewhole space into a list of objects with corresponding
properties, such as location and uncertainties (Thrun etal. 2005). For example, a 2D space can be generalized
into an “occupancy grid” (Elfe 1989). The occupancy
grid represents a dense 3D map with a finite array ofpoints. Any additional measurements or data sources
over the same space can then relate to one or several

Page 3 of 13

grid points. A similar concept of sparse discretized mapping has also been applied into our approach to model
topographic features. A 3D surface in space can be approximated using a finite set of known points, on which
any given point can be predicted via a wide variety ofinterpolation techniques. Some of these techniques arecommonly used in geosciences, including Kriging
models (Krige 1951) and other types of regression
methods (Williams 1998; Dumitru et al. 2013).
There is also a realization that it is equally important,
if not more so, to accurately represent the uncertainty
function across the whole space. Spatial uncertainty
modeling is a key element in the regression-based prediction processes. Moreover, the spatial uncertainty
functions may change over time due to the underlying
geomorphic processes, and this can be evaluated with astochastic estimator, such as a Kalman filter (Kalman
1960). For example, Mardia et al. (1998) and Cortes(2009) suggested the combination of spatial and temporal modeling by using Kriging and Kalman filters.
Here, our focus is on the spatio-temporal uncertainty
model itself, instead of any specific interpolation
method. Our approach produces an efficient, accurate,
and robust uncertainty model that opens the door to theintegration of legacy data and new sensors, and provides
more definitive measures of landform and landscapeevolution from a variety of sources. A major benefit andadvance of our adopted approach is that an optimum
interpolation method can be applied to this model,
which would estimate or predict the elevation and theassociated uncertainty for any location in the entire region, at any given time in history or even in the future.

Methods/Experimental
Spatio-temporal uncertainty model

Here, we consider topographic uncertainty as a stochastic process that takes into consideration spatial and temporal variations. The spatial uncertainty can be modeled
as a Gaussian process, which tends to vary across the region of interest, and yet often has local correlations. A
higher level of spatial correlation may be expected fromsmooth surfaces, whereas sudden elevation changes,
such as a steep cliff, will result in lower correlation. A
covariance function can describe this type of uncertainty
model and is required by most interpolation or prediction techniques. Although it could be challenging to obtain the complete covariance function, especially if suchfunction is also changing over time, following Williams(1998), we estimate a covariance matrix over a finite setof hypothetical anchor points. The covariance matrix is
updated over time, based on input from a variety of sensors, such as laser scanners and global navigation satellite system (GNSS) survey. Although the filter does

Wasklewicz et al. Progress in Earth and Planetary Science (2017) 4:32

benefit from direct observations whenever available, sensors are not expected to make direct observation of theanchor points. For example, a terrestrial laser scanner(TLS) provides a high-resolution scan of the region, but
the point cloud is not guaranteed to overlap the hypothetical anchor points. Static GNSS surveys in historical
data, on the other hand, may only be available on a few
points in a region, which are even less likely to overlap
with the anchor points. The lack of coincident data withthe anchor does not preclude the use of this method, rather this method does not require any direct sensor observations in our experimental design or in any future
field applications of this methodology.
Anchor point distributions are instead dependent
upon the complexity of the topographical features,
which is independent of any specific sensor measurements. However, sensor observation made at any given
location within the region can still be used to update theuncertainty model of the neighboring anchor point(s).
To achieve that, we further assume that more points will
be used to model an area with lower spatial correlation,such that the elevation at a given location can be interpolated or extrapolated with neighboring anchors with
sufficient accuracy. Thus, a surficial feature with topographic complexity is covered with densely populated
anchor points, whereas less topographically complex
areas are covered with a sparse scatter of anchor points.
The uncertainty model and associated anchor points
not only account for the spatial uncertainty, but also
consider the function of time. A stochastic estimator
traces the evolution of such randomness, in which theanchor points form a space of “states”. Thus, the elevation and uncertainty of any given point at any time canbe approximated with a combination of these states. In asoftware simulation, we created a landform of 150 by
150 m, and simulated the elevation change over thecourse of 30 years. Both the landform and the elevation
change have non-linear surfaces (Fig. 1). To demonstrate
the concept of using anchor points, we first formulated

Page 4 of 13a 5 by 5 array, which has sufficient density for a discretized representation of the simulated landform. However,anchor points do not necessarily form any regular shapein general.
Data from a TLS, an ALS, a 2D aerial camera, a stereo
vision system, structure from motion SfM-MVS from alimited number of historic aerial photographs takenfrom a nadir-looking camera, static and Kinematic
GNSS surveys are considered in our simulations(Table 1). These sensors differ greatly in accuracy andresolution. We assume all the sensor errors can be projected onto the vertical direction, and thus only the vertical errors are of concern in developing an objective
comparison among the various sensors.
Gaussian error models are used for GNSS positioning
in this simulation. Naturally, the static surveys are moreaccurate. However, the Kinematic GNSS survey will
gather more data points (Table 1).
The uncertainty in laser ranging measurements can bemodeled with two main components, angular and radial.
Although the complete ranging error model may besophisticated (Thrun et al. 2005), it is often sufficient toassume that in practice, both error sources follow anormal distribution. The underlying mechanisms are independent from each other (Glennie 2007). Figure 2 illustrates a typical application of laser scanner, whichapplies to both TLS and ALS. The laser scanner is elevated from the ground by h, and it scans a location on
the ground at a line of sight distance r. Errors in thehorizontal (azimuth) and vertical (elevation) angles result from laser beam width and the precision in orientation measurements and follow a 2D normal distribution(red ellipse in Fig. 2). Since the focus of this work is noton the high-fidelity simulation of ALS data, the comprehensive error model (Glennie 2007) is abstracted andrepresented with two components: a radial error σr, anda vertical angular error σθ. A target observed at distancer has the total uncertainty of σv. Compared to an ALS, a
TLS is closer to the ground, produces a denser point

Fig. 1 Simulated landform (left) and accumulated changes over 30 years (right). The magnitude of change is approximately ± 5 m

Wasklewicz et al. Progress in Earth and Planetary Science (2017) 4:32

Page 5 of 13

Table 1 Sensor resolution and error
Sensor name

TLS

ALS

2D camera

SfM

Static GNSS

Kinematic GNSS

Vertical error standard deviation in simulated data (m)

0.02

0.1

N/A

3.3

0.04

0.08

Sample Density (pt/km2)

1.1 × 107

7 × 103

1.1 × 105

1.1 × 105278

2.2 × 104

The values presented herein are representative of our instrumentation and our experience to represent sensor quality. The large vertical error is observed with
SfM since the displacement between camera perspectives is short in the simulation and is represented of legacy aerial photography. Sample density represents
the density of measurements used to update the estimator in this simulation. It does not necessarily represent the density of raw data available from thesesensors (such as TLS and ALS)

cloud, and has more accurate ranging measurements.
The ALS, on the other hand, scans the ground from ahigh altitude, in which case the vertical error is less sensitive to the angular uncertainty.
Data from a 2D camera image does not provide direct
observation of elevation, as can be seen in Table 1. However, 2D images can be related to the 3D landform via a
3D-2D projection model (Hartley and Zisserman 2003).
Let vector [u, v]T represent the 2D location of a point L
on the image of a landform, observed using a camera
Twith lens focal length f, and vector xCL ; yCL ; zCL
be the3D location of point L, as observed in the camera frame(C frame). The 3D-2D relationship can be modeled with
" #
 
f xCL
u
¼ C C ;v
z L yL
which is implemented in the estimator. Any change inxCL and yCL are thus observable in 2D images. The landform elevation change is observed on the z axis in a Global frame G, ΔzG
L . As shown in Fig. 3, point L is locatedat X G,
anditisbeing observed by a camera located at
L
XG
,
both
in
the
globalframe. The elevation error ΔzG
L isC
related to the observation error made in the camera
frame ΔxCL . ΔxCL can be geometrically decomposed into
errors in the horizontal direction, ΔxG
L and the verticaldirection, ΔzG
intheglobalframe.
The
camera associLated with the aerial photography is typically locatedabove the landform (Fig. 3). The vertical observable ΔzG
L
is thus smaller than the horizontal component ΔxG
L inthis geometry. Therefore, it tends to be less effective(Table 1) when elevation changes are detected with an
overhead image collected with an airborne camera.
The aforementioned projection model should not be
confused with SfM-MVS, which operates on the same
principle of multi-view geometry (Hartley and Zisserman2003), as does stereo vision or triangulation. A camera
in motion (with a limited overlapping nadir-looking images to simulate a historic aerial photogrammetric approach) will provide repeated observations of the same3D landform [u, v]T1. . n, collected from n perspectives(unstructured photographs). In this process, the camera
orientation and location are known at all n perspectives,
referenced to frame G. The 2D images, [u, v]T1. . n willthen be converted into unit vectors in frame G, UG
1. . n.
At the ith perspective, UiG points from the camera posGG
ition X G
C;i to the landform position X L . Since X L does
not change over the different perspectives, it is solved byGusing UG
1. . n and X C;1::n via an optimization process(Hartley and Zisserman 2003). Figure 4 illustrates thegeometric relationship between the landform and thecamera on perspectives 1 and n. The uncertainty of X G
Lwill also be dependent on the accuracy of UG
1. . n andG
XG
C;1::n , and the geometric relationship between X C;1::nand X G
L . We conservatively assume a short displacement
between camera perspectives in this experiment, which
results in greater vertical errors in SfM.
The stochastic estimator uses measurements that become available sequentially over time. Let a state vector
x represent the vertical errors on the 25 anchors, locatedat mx[1..25], and a 25 by 25 matrix P for the covarianceamong states. x and P are initialized with someFig. 2 Laser scanner error model. θ: vertical angle, r: range, σr: standard deviation of radial ranging error, σθ: standard deviation of vertical angularerror, σv: standard deviation of total vertical error

Wasklewicz et al. Progress in Earth and Planetary Science (2017) 4:32

Page 6 of 13

Fig. 3 2D camera image error model in a typical aerial photography application

measurements, such as a TLS scan at time tk-1, and thenpropagated to another time tk using a Brownian process.
xk;k−1 ¼ Fxk−1;k−1 þ μ
P k;k−1 ¼ FP k−1;k−1 F þ Q;
where F is a matrix describing the dynamic relationshipbetween states, which will be the identity matrix in thiscase. μ is the process noise component, corresponding
to the Brownian process, of which the covariance matrix
is defined with Q.
Without additional measurements, the matrix P, willgrow as a function of time. The Brownian process may
be established with the best-known model of terrainchange, and will only be used to predict the increase of
uncertainty if the model is valid for the underlying geomorphic process.

The fundamental concept can be illustrated using justtwo of the anchor points, located at mx[1] and mx[2], respectively. The initial uncertainty of both points can be
represented with two circles at t0 (Fig. 5). As time propagates to t1, the uncertainty grows for both points, illustrated with greater circles. A new measurement
becomes available at t1, made at location mz[t1]. The
uncertainty of this measurement is represented with anellipse. This measurement is then used to update bothanchor points. Right after this update, at t1+, the uncertainty of both points is reduced to ellipses. Since theupdate is closer to mx[2], the corresponding ellipse becomes smaller than that of mx[1]. Without additional
measurements, both ellipses are propagated to t2, which
result in greater ellipses. A second measurement wasmade at t2, which occurred much closer to mx[1] thistime. It effectively reduced the size of the ellipses, more
so on mx[1] than mx[2].

Fig. 4 SfM error model with two perspectives. A long displacement example is shown on the left, with the error of landmark location illustrated
in a purple ellipse. On the right is a short displacement example, where the error is greater in the vertical direction

Wasklewicz et al. Progress in Earth and Planetary Science (2017) 4:32

Page 7 of 13

xk;k ¼ xk;k−1 þ Ky
P k;k ¼ ðI−KH ÞP k;k−1

Fig. 5 Two anchor points located mx[1] (blue) and mx[2] (red) showwhat happens when new measurements become available overtime. Uncertainty in both points are represented by circles and thecircles increase and decrease in size with increasing or decreasing
uncertainty, respectively

The geometric relationship between anchor points andany measurement can still be represented or approximated with a linear/linearized model, which can beestablished by using an interpolation method. Noticethat the interpolation model used in this step estimatesthe covariance function, instead of predicting the landform as does Kriging or polynomial regression. The bicubic interpolation method introduced in Keys (1981) has
been widely used in the computer vision community, andit has been adopted in digital elevation surface modeling(Dumitru et al. 2013). With bicubic interpolation, a linearized relationship between a given point and the anchor
points, denoted with matrix H, can be easily estimated.
This approach handles non-linear landforms, and remains
relatively computationally efficient.
The elevation at mz now has a direct measurement z,
and a prediction from the estimator states, Hxk; k − 1.
Any disagreement between them is likely a result of uncertainties associated with both and new informationthat can be observed as “innovation” on the states,
y ¼ z−Hxk;k−1 ;with a de-facto standard estimation method, the Kalmanfilter (Kalman 1960), x and P are updated with a gain K:−1
K ¼ P k;k−1 H ′ HP k;k−1 H ′ þ Rwhich allows us to update the states x and covariance P
using this new information:

After this step, the states are further propagated onto
time tk + 1, with the integration of another measurement(such as aerial photo). The time history of states x andcovariance matrix P are the outcome of this approach.
Notice that only a few key steps of the estimator are
highlighted here for the sake of conciseness. For details
on the derivation and implementation, please refer to(Kalman 1960) and (Smith et al. 1962).
At any given point in time, x and P are used as input
to predict the elevation and uncertainty of any given
point on the landform with a standard interpolation
method. An accurate and consistent covariance matrix P
is the key to an optimum interpolation process.
Consistency of covariance is defined based on how well
P describes and over-bounds the actual error in states x.
To verify the estimation, x will be compared against a
truth reference xref used in the software simulation.
With field data, the truth reference can be generated
from a high-resolution, high-fidelity sensor such as TLS.
The variance and covariance of error Δx = x − xref, is expected to be closely bounded by P at each time step tk.
Since it may be difficult to visualize the comparison between matrices, the standard error σΔx[k] will be compared against an uncertainty level for this step, by using
the root mean square of all the diagonal elements in P.
qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ﬃ
σ est ½k  ¼ diag P k;k
In the software simulation, the actual elevation change
is not provided to the estimator. Rather, an overbounding Brownian process is used, and we expect tohave σ est ≥σ Δx in the history of 30 years.
Physical modeling experiment

A physical model is developed in the Terrain Analysis
Lab at East Carolina University as a further means ofvalidating our uncertainty measures. The controlled
simulation is designed around a fan-shaped surficial features (i.e. an alluvial fan or washover fan) built in astainless-steel stream table (3.7 m long by 1.8 m wide)using coarse sand on a white sheet of plotter paper toreduce reflectivity of the stainless-steel table (Fig. 6).
The fan-shaped feature is 0.974 m wide and 0.392 m inwith a relief of 0.05 m from the bottom of the stream
table (Fig. 6).
A Leica P20 laser scanner was inverted and mounted
on an aluminum swing set (Lisenby et al. 2014). Targets
mounted on the wall of the lab and within the stainlesssteel table as well as the corner points of the table wereused to register the TLS data (Fig. 6). CartesianWasklewicz et al. Progress in Earth and Planetary Science (2017) 4:32

Page 8 of 13

Fig. 6 The experimental fan model. a Location of control points and cameras. b SfM image of the fan on the plotter paper with the circle targets(image rotated 90 CCW to A). c A topographic map of the experimental fan. d Elevation difference from t1 to t0

coordinates associated with TLS point clouds provide locations for the control points used in the SfM. A single
scan from the nadir looking position of the scanner captured the entire simulated fan surface in t0 and t1. SfM
data were collected with a Ricoh GR II digital camera.
The camera was positioned at two different heights at 32
locations looking obliquely and inward at the fan surface(Fig. 6). A total of 64 images were used to generate apoint cloud with the aid of Agisoft Photoscan software.
A combination of 14 12-bit targets from the Agisoft software were printed on sticky back paper and adhered tothe stainless-steel table, along with 10 small circle targets with known dimensions. Targets and circles wereused to aid in the photo alignment process.
The original landform (t0) was modified (t1) by adding
approximately 2 cm of coarse sand onto a segment ofWasklewicz et al. Progress in Earth and Planetary Science (2017) 4:32

the fan (Figs. 6 and 7) to emulate fan segment aggradation over an arbitrary time-period (for example, t1 is set
to 30 years). TLS and SfM-MVS data are collected atboth t0 and t1. TLS at t1 is used as a truth reference. We
initialize the estimator at t0 with TLS and SfM-MVS,
propagate and update at t1 using SfM-MVS. The estimated σ est ½1 is compared against the actual error σΔx[1]in this case.

Results and discussion
The estimator in the software simulation is initializedwith a TLS scan. Subsequently, five additional sensor
data sets become available over 30 years, to observe thegeomorphic change. A general pattern evolves whereby
the estimated uncertainty provides an over bound of theactual error (Fig. 8), i.e., σ est ≥σ Δx . Both σ est and σΔxgrow over time during the simulation, and are dramatically reduced by some measurements, such as ALS, SfMMVS and Kinematic GNSS surveys. Other sensors, suchas 2D photos and static GNSS, are less effective, as expected. Although the SfM-MVS measurements are notnearly as accurate as static GNSS survey, it offers amuch higher point density and therefore, is also able toreduce the error more effectively. The over-bounding associated with the static GNSS has more to do with thelack of data produced with this technique than the accuracy of the data collected. The truth reference used inthis software simulation can be used to verify the estimation in x (Fig. 9).
The estimator is also applied to the changes associatedwith the lab-simulated fan model. As this terrain is moresophisticated, a 20 by 20 mesh of anchors is used. The
actual change xref is obtained by comparing the two TLS
scans (Fig. 10). The estimated change is observed on the

Page 9 of 13

anchor points, by using a less accurate and lowerresolution SfM-MVS update. At t1 = 30 years, σ est ¼ 5:6
mm and σΔx = 3.7mm, σ est is greater than, and yet reasonably close to σΔx. The difference is below the detectable limits for topographic change associated with thedata source instrumentation (Wester et al. 2014, Staley
et al. 2014; Wasklewicz and Scheinert 2015).
The software simulation and physical modeling experiment both show that the estimator can provide a reliable
and consistent measure of uncertainty in legacy data. Although some sensors are more desired than others, theavailability of sensors in legacy data and sometimes also
in new data collection campaigns dictates that data fromall sensors must be incorporated. The estimator willallow us to retrieve information from sensors of low accuracy, low observability (2D images) or low spatial
resolution (static surveys). We further emphasize thatthe discussed estimator is not limited to just the Kalmanfilter. Other Bayesian estimators, such as particle filters(Montemerlo and Thrun 2003), are applicable as well.
Similarly, bicubic interpolation is not the only choice forinterpolation either. Other methods that can yield a linearized relationship between measurements and states arealso valid options, which will provide various types ofadvantages (Dumitru et al. 2013).
The statistical consistency of uncertainties in all the
sensor measurements and estimates are critical in ourapplication. In practice, repeated measurements could
also lead to a covariance matrix P that is too small(Julier and Uhlmann 2001), such as high-density datafrom SfM-MVS. As x and P can be used to predict theelevation and uncertainty of the whole landform, anoverly optimistic P can easily cause unexpected errors in
the prediction. For instance, few erroneous points in x

Fig. 7 Elevation model from TLS point cloud of the simulated fan model at t1

Wasklewicz et al. Progress in Earth and Planetary Science (2017) 4:32

Page 10 of 13

Fig. 8 Uncertainty of software simulated landform over 30 years, estimated sigma σ est vs. actual error σΔxwith small variances would dominate the prediction andgenerate a bias in a neighborhood. Overly optimistic P
also would prevent x from being updated with othersensors. Additional steps will be taken in future work toguarantee the consistency of P.
The successful outcomes and corroboration in our initial software and physical modeling simulations, the ability
to adopt and use multiple estimators and interpolation
techniques, and an ability to be employed across a varietyof environmental settings provides a new alternative to integrating legacy and HRT data.
Applications of existing and planned methodology tofield research will help extend our understanding of the

dynamic landform and landscape changes that take placethrough space and time. At present, earth science researchers conducting studies with the aid of HRT data
have placed a great deal of emphasis on methods andtechniques that provide rapid, accurate, and spatiallycontinuous topographic data. HRT provide the user with
unparalleled information of short-term changes in processes and forms. Subtle topographic changes that earthscientists must measure in the future, for example those
associated with climate change (Lane 2013), require detailed data that can be precisely acquired and accurately
assessed to address scientific questions and inform policy makers. However, the sole use of HRT data limits theFig. 9 Actual vs. estimated terrain change in year 7 (left), year 19 (middle) and year 31 (right). Upper: actual terrain change xref; lower: estimationin states x

Wasklewicz et al. Progress in Earth and Planetary Science (2017) 4:32

Fig. 10 Actual vs. estimated terrain change in lab-simulated fanmodel. Upper: actual terrain change xref; lower: estimation in states x

temporal framework for which earth scientists can consider these changes because these data sets have only
been recently developed. Legacy topographic data allow
us to extend this temporal perspective of change backward in time to better inform us of variations on decadal
time-scales and in some instances, as far back as a halfcentury or more (Corsini et al. 2009; Carley et al. 2012;
James et al. 2012; Schaffrath et al. 2015). While we appreciate that this represents a miniscule amount of timefrom a geological perspective, the ability to use legacy
topographic data increases the temporal sample size
earth scientists are working from and allows us to makemore authoritative statements about landform evolution
than the shorter-term view supplied by research solely
using HRT data.
Scientists need to be able to definitively measure theenvironmental changes and communicate these measures to a broad audience of possible users that include
their peers, public officials, and the community. Expressing the uncertainty associated with these measurementsor reducing the uncertainty in the data production
process is a critical first step to achieving more definitive
measures and communication. Here, our uncertainty
model, associated anchor points, and stochastic estimator clearly highlight that spatial uncertainty can be
accounted for across the time-scales of the legacy and
HRT data. Each elevation point/grid has an associatedPage 11 of 13

error, which can be used to produce better maps wheninterpolating to raster data and provide better measurement tolerances as they relate to topographic changes
and sediment transport.
The ability to accurately assess spatially variable uncertainty across a broad range of temporal scales also provides a means to exchange information between
researchers, engineers, policy makers, and environmental
managers and planners. The longer-term understandingof topographic change and sediment transport and theuncertainties associated with these dynamics permits
greater understanding of issues such as: water quality,
habitat loss, and risks to built-environments and human
lives from environmental hazards. Furthermore, an ability
to measure the level of data uncertainty associated withpast landform and landscape dynamics will help advancelandscape evolution models because of the ability to corroborate existing landscape evolutionary models with thelonger-term record of topographic change derived fromlegacy topographic data sets and recent HRT data. This isa necessary advancement as we may come to rely on landscape evolution models as an integral tool in the decisionand management schemes as we face climate change andincreased scenarios of extreme weather.

Conclusions
Our software simulation and physical modeling experiments provide a new approach to measuring topographic data uncertainty where legacy data from avariety of data sources can be integrated with HRT data
to expand the time-scales of topographic change detection. As anticipated, the difference in the actual and expected errors of our HRT physical model experimentwas quite small (< 2 mm). Current instrumentation andfield methods often have a higher minimal level of detection, so this value is quite acceptable for HRT data
sources and represents a value comparable to or lower
than most uncertainty measures found in current research
exploring topographic change detection. The uncertainty
model, associated anchor points, and stochastic estimator
was further applied to a software simulation whereby avariety of remote-sensing data sources were used to simulate data capture from legacy data sources. Our findings
show the estimated error coincides with the actual error
using certain sensors (Kinematic GNSS, ALS, TLS, and
SfM-MVS). Data from 2D imagery and static GNSS did
not perform as well at the time the sensor is integrated
into estimator. Nevertheless, the software simulationshows the approach can be used to estimate the error associated with all elevation values in the legacy and HRT
data over the time-period of the simulation.
Our findings show a strong potential for this technique to be applied in a variety of field settings. We anticipate further development of the uncertainty model inWasklewicz et al. Progress in Earth and Planetary Science (2017) 4:32

future research. Our goal is to expand the capability todensify the anchor points in areas that are more topographically complex and test by overfitting the surface
how well the uncertainty model is performing in theselocations. At present, it is not clear if the use of sparse
points (for example, techniques that rely on interpolationfrom topographic maps) under sample the topographic
surface and therefore, do not accurately represent the surface or the spatially variable uncertainty in the measurement of topographic change. We intend to expand thecapabilities of the model by conducting further software
simulations that analyze and visualize topographic changeover time using the varied densification of anchor points.
The advances in the uncertainty model will then be further evaluated using real-world data from a variety ofsources. For example, ALS data can be used in a fashion
to simulate different types of data by varying the data
density, data noise, etc. to examine known differences inthe data and the capability to measures these data perturbations using the refine uncertainty model. Once these
items have been worked out our intent is to apply the uncertainty model to the real-world topographic and environmental changes within various environmental settings
that include more topographically complex areas likemountainous environments versus less topographic complex locales such as barrier islands.
Our existing uncertainty model shows clear evidence
that the temporal component of spatially variable uncertainty of legacy topographic data sets can be measured
in both of our simulations. A capacity to integrate legacy
topographic data expands earth scientist’s ability tounderstand topographic changes and sediment transport
over longer time-scales. This should lead to more definitive measurements and answers regarding how environments are changing over the past 50 to 100 years(dependent upon the temporal availability of topographic
data sources). More conclusive measures and findings
stem from our ability to assess uncertainty across the
spectrum of legacy data in relation to more recent HRT
data using techniques like we have documented in thisresearch. New results from these approaches should provide valuable information to the broader scientific community and society. Societal benefits will likely berecognized in the form of delivering relevant information that show the potential range for environmental
changes over time-scales demanded by managers andpolicy makers.
Abbreviations
ALS: Airborne laser scanning; DEM: Digital elevation model; DEM-of-difference;
GNSS: Global navigation satellite system; HRT: High-resolution topography; SfMMVS: Structure from Motion-Multi-View Stereo; SIE: Survey and interpolation
error; TIN: Triangular irregular networks; TLS: Terrestrial laser scanning
Acknowledgements
We thank Kailey Adams and Robert Rankin for their assistance in our experiments.

Page 12 of 13

Funding
This work was funded by an East Carolina University Interdisciplinary Research
Collaboration Award entitled “Development of Scalable Sensor Payload for Rapid
High Resolution Topographic Scanning via Micro Aerial Vehicles”. Funding for thisaward was supplied by the Division of Research and Graduate Studies inpartnership with the Divisions of Academic Affairs and Health Sciences.
The Terrain Analysis Lab and ECU Department of Geography, Planning, andEnvironment, provided the physical modeling facilities and equipment tosimulate field measurements in a controlled setting.
Authors’ contributions
TW proposed the topic, helped conceive the study design, analyzed thedata, helped with the interpretation of the data, co-wrote the paper, and iscorresponding author. ZZ helped conceive the study design, carried out theexperimental study, analyzed the data, interpreted the data, designed figures,and co-wrote the paper. PG helped in the interpretation of the data andco-wrote the paper. All authors read and approved the final manuscript.
Competing interests
TW, ZZ, and PG declare that they have no competing interest.

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.
Author details
1
Department of Geography, Planning, and Environment, East Carolina
University, Brewster Building A227, Greenville, NC 27858, USA. 2College ofEngineering and Technology, East Carolina University, Science andTechnology Complex Suite 100, Greenville, NC 27858, USA.
Received: 28 January 2017 Accepted: 28 September 2017

References
Brasington J, Langham J, Rumsby B (2003) Methodological sensitivity of morphometricestimates of coarse fluvial sediment transport. Geomorphology 53:299–316
Brasington J, Rumsby BT, McVey RA (2000) Monitoring and modellingmorphological change in a braided gravel-bed river using high resolution
GPS-based survey. Earth Surface Processes and Landforms 25: 973-990
Carley JK, Pasternack GB, Wyrick JR, Barker JR, Bratovich PM, Massa DA, Reedy GD,
Johnson TR (2012) Significant decadal channel change 58–67 years post-dam
accounting for uncertainty in topographic change detection between contour
maps and point cloud models. Geomorphology 179:71–88
Corsini A, Borgatti L, Cervi F, Dahne A, Ronchetti F, Sterzai P (2009) Estimating
mass-wasting processes in active earth slides–earth flows with time-series ofhigh-resolution DEMs from photogrammetry and airborne LiDAR. Nat
Hazards Earth Syst Sci 9:433–439
Cortes J (2009) Distributed kriged Kalman filter for spatial estimation. IEEE Trans
Autom Control 54:2816–2827
Crowell M, Leatherman SP, Buckley MK (1991) Historical shoreline change: error
analysis and mapping accuracy. Journal of coastal research, 839-852
Dumitru PD, Marin P, Dragos B (2013) Comparative study regarding the methodsof interpolation. In: Recent advances in geodesy and Geomatics engineering,
pp 45–52
Elfe A (1989) Using occupancy grids for mobile robot perception and navigation.
Computer 22:46–57
Fonstad MA, Dietrich JT, Courville BC, Jensen JL, Carbonneau PE (2013) Topographic
structure from motion: a new development in photogrammetric measurement.
Earth Surf Process Landf 38:421–430
Glennie C (2007) Rigorous 3D error analysis of kinematic scanning LIDAR systems.
J Appl Geodesy 1:147–157
Glennie CL, Hinojosa-Corona A, Nissen E, Kusari A, Oskin ME, Arrowsmith JR,
Borsa A (2014) Optimization of legacy lidar data sets for measuring near-field
earthquake displacements. Geophys Res Lett 41:3494–3501
Gomez C, Hayakawa Y, Obanawa H (2015) A study of Japanese landscapes using
structure from motion derived DSMs and DEMs based on historical aerial
photographs: new opportunities for vegetation monitoring and diachronicgeomorphology. Geomorphology 242:11–20
Hartley R, Zisserman A (2003) Multiple view geometry in computer vision.
Cambridge University Press, CambridgeWasklewicz et al. Progress in Earth and Planetary Science (2017) 4:32

Heritage GL, Milan DJ, Large AR, Fuller IC (2009) Influence of survey strategy andinterpolation model on DEM quality. Geomorphology 112:334–344
James LA, Hodgson ME, Ghoshal S, Latiolais MM (2012) Geomorphic change
detection using historic maps and DEM differencing: the temporal dimension of
geospatial analysis. Geomorphology 137:181–198
Julier SJ, Uhlmann JK (2001) General decentralized data fusion with covarianceintersection (CI). In: Handbook of data fusion. CRC Press, Boca Raton
Kalman RE (1960) A new approach to linear filtering and prediction problems.
J Basic Eng 82:35–45
Keys R (1981) Cubic convolution interpolation for digital image processing. IEEE
Trans Acoust Speech Signal Process 29:1153–1160
Krige DG (1951) A statistical approach to some mine valuations and allied problems
at the witwatersrand. Master’s thesis. University of Witwatersrand, WitwatersrandLane SN (2013) 21st century climate change: where has all the geomorphology
gone? Earth Surf Process Landf 38:106–110
Lane SN, Westaway RM, Murray-Hicks D (2003) Estimation of erosion anddeposition volumes in a large, gravel-bed, braided river using synoptic
remote sensing. Earth Surf Process Landf 28:249–271
Lisenby PE, Slattery MC, Wasklewicz TA (2014) Morphological organization of asteep, tropical headwater stream: the aspect of channel bifurcation.
Geomorphology 214:245–260
Mardia KV, Goodall C, Redfern EJ, Alonso FL (1998) The Kriged Kalman filter. TEST
7:217–285
Micheletti N, Lane SN, Chandler JH (2015) Application of archival aerial
photogrammetry to quantify climate forcing of alpine landscapes.
Photogramm Rec 30:143–165
Milan DJ, Heritage GL, Hetherington D (2007) Application of a 3D laser scanner in
the assessment of erosion and deposition volumes and channel change in a
proglacial river. Earth Surface Processes and Landforms 32:1657-1674
Milan DJ, Heritage GL, Large AR, Fuller IC (2011) Filtering spatial error from
DEMs: implications for morphological change estimation. Geomorphology,
125: 160-171
Montemerlo M, Thrun S (2003) Simultaneous localization and mapping withunknown data association using FastSLAM. Proceedings of IEEE International
Conference on Robotics and Automation, Taipei
Oskin ME, Arrowsmith JR, Corona AH, Elliott AJ, Fletcher JM, Fielding EJ, Gold PO,
Garcia JJG, Hudnut KW, Liu-Zeng J, Teran OJ (2012) Near-field deformationfrom the el mayor-Cucapah earthquake revealed by differential LIDAR.
Science 335:702–705
Pelletier JD, Brad Murray A, Pierce JL, Bierman PR, Breshears DD, Crosby BT, Ellis M,
Foufoula-Georgiou E, Heimsath AM, Houser C., Lancaster, N (2015) Forecasting
the response of Earth's surface to future climatic and land use changes: A
review of methods and research needs. Earth's Future, 3(7)220-251
Schaffrath KR, Belmont P, Wheaton JM (2015) Landscape-scale geomorphic
change detection: quantifying spatially variable uncertainty andcircumventing legacy data issues. Geomorphology 250:334–348
Smith GL, Schmidt SF, McGee LA (1962) Application of statistical filter theory to
the optimal estimation of position and velocity on board a circumlunar
vehicle. National Aeronautics and Space Administration, Washington, D.C
Spencer T, Naylor L, Lane S, Darby S, Macklin M, Magilligan F, Möller I (2017)
Stormy geomorphology: an introduction to the special issue. Earth Surf
Process Landf 42:238–241
Staley DM, Wasklewicz TA, Kean JW (2014) Characterizing the primary material
sources and dominant erosional processes for post-fire debris-flow initiation
in a headwater basin using multi-temporal terrestrial laser scanning data.
Geomorphology 214:324–338
Tarolli P, Sofia G (2016) Human topographic signatures and derived geomorphic
processes across landscapes. Geomorphology 255:140–161
Thrun S, Burgard W, Fox D (2005) Probabilistic robotics. MIT press, Boston
Wasklewicz T, Scheinert C (2015) Development and maintenance of atelescoping debris flow fan in response to human-induced fan surface
channelization, Chalk Creek valley natural debris flow laboratory, Colorado,USA. Geomorphology 252:51–65
Wasklewicz T, Staley DM, Reavis K, Oguchi T (2013) Digital terrain modeling. In:
Shroder J, Bishop MP (eds) Treatise on geomorphology, vol 3. Academic
Press, San Diego, pp 130–161
Wasklewicz TA, Hattanji T (2009) High-resolution analysis of debris flow–Induced
Channel changes in a headwater stream, Ashio Mountains, Japan. Prof Geogr
61:231–249

Page 13 of 13

Wester T, Wasklewicz T, Staley D (2014) Functional and structural connectivity
within a recently burned drainage basin. Geomorphology 206:362–373
Wheaton JM, Brasington J, Darby SE, Sear DA (2010) Accounting for uncertainty
in DEMs from repeat topographic surveys: improved sediment budgets.
Earth Surf Process Landf 35:136–156
Williams CKI (1998) Prediction with Gaussian processes: from linear regression tolinear prediction and beyond. Learn Graphical Models 89:599–621