Next Article in Journal
Use of SSU/MSU Satellite Observations to Validate Upper Atmospheric Temperature Trends in CMIP5 Simulations
Previous Article in Journal
User Validation of VIIRS Satellite Imagery
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Technical Note

Estimating Horizontal Displacement between DEMs by Means of Particle Image Velocimetry Techniques

Department of Architectural Graphic Expression and Engineering, University of Granada, 18071 Granada, Spain
*
Author to whom correspondence should be addressed.
Remote Sens. 2016, 8(1), 14; https://doi.org/10.3390/rs8010014
Submission received: 31 July 2015 / Revised: 12 December 2015 / Accepted: 16 December 2015 / Published: 24 December 2015

Abstract

:
To date, digital terrain model (DTM) accuracy has been studied almost exclusively by computing its height variable. However, the largely ignored horizontal component bears a great influence on the positional accuracy of certain linear features, e.g., in hydrological features. In an effort to fill this gap, we propose a means of measurement different from the geomatic approach, involving fluid mechanics (water and air flows) or aerodynamics. The particle image velocimetry (PIV) algorithm is proposed as an estimator of horizontal differences between digital elevation models (DEM) in grid format. After applying a scale factor to the displacement estimated by the PIV algorithm, the mean error predicted is around one-seventh of the cell size of the DEM with the greatest spatial resolution, and around one-nineteenth of the cell size of the DEM with the least spatial resolution. Our methodology allows all kinds of DTMs to be compared once they are transformed into DEM format, while also allowing comparison of data from diverse capture methods, i.e., LiDAR versus photogrammetric data sources.

Graphical Abstract

1. Introduction

A lot of geomatic applications use DTMs as the basis to derive other products such as cartography, hydrological features, engineering solutions, urban planning, environmental studies, time series analyses, reports, statistics, etc. The precision and accuracy of DTMs, required for deriving the above-mentioned products, however, may vary substantially. The data capture methodology is one key to explain the precision and accuracy of DTMs and how they influence the derived product. There is a need for reliable measures of accuracy allowing comparative evaluations of DTMs from different sources, e.g., LiDAR versus photogrammetry. Most studies thus far focus on the height variable [1,2]. As the easiest variable to statistically compute and interpret, such studies focus on vertical accuracy. Some authors [3,4] studied the DEM accuracy as a function of some features (topographic index, hydrological objects). There are no reports to date of horizontal accuracy of DTMs or DEMs linked with the derived features, most likely because any such measure, if it exists, would hardly be applicable to other features [5]. More recently, algorithms were defined to link the differences between horizontal accuracy and the derived parameters of features; such algorithms could be based on a triangulated irregular network (TIN) [6,7] or DEM grid models [8]. The first approach generalizes the separation between homologous contours as the measure of horizontal discrepancy for the whole DEM, whereas the second compares several cells in two or more DEMs to estimate the horizontal discrepancy.
We can define the vertical accuracy for a DEM (grid structure composed of rows and columns) as the magnitude indicating how close the height of a point is, located by its horizontal position (row and column), on a DEM with respect its true value. If we have a reference DEM (RefDEM) as the ground truth we could estimate the vertical accuracy from any other DEM (CompDEM) just by computing the differences between the corresponding cell’s height from both DEMs. The vertical accuracy can be seen as the matrix: accMat = RefDEM − ComDEM. Horizontal accuracy is a little bit more complicated to understand because it may be difficult, given a point (RefP) on a RefDEM, to find its homologous point (HomP) on a CompDEM, except, for example, when RefP and HomP are a building corner, but that is not usually the case. We could understand the horizontal accuracy of a CompDEM with respect to its RefDEM; imagining the case where the CompDEM is identical to the RefDEM, then their contours would be the same. However, if we move the identical CompDEM horizontally, then the homologous contours from the CompDEM and RefDEM would not match and the discrepancy between the homologous contours indicates the horizontal discrepancy, which, consequently, is a measure of the horizontal accuracy. It is possible for each cell in the CompDEM to be horizontally displaced in any direction, so that homologous contours cross each other as shown in Figure 1, where the area between homologous contours has been filled with gray color.
Figure 1. Homologous contours (coming from homologous DEM: RefDEM and CompDEM) crossing each other and filling the area with gray color. The black color is a river piece.
Figure 1. Homologous contours (coming from homologous DEM: RefDEM and CompDEM) crossing each other and filling the area with gray color. The black color is a river piece.
Remotesensing 08 00014 g001
Both vertical and horizontal accuracy are important features to keep in every DEM because error in vertical and horizontal positioning could have bad consequences when, for example, computing flood risk areas. If we miscompute flood areas as a consequence of inaccurate DEMs, it would be possible to authorize, by error, the construction of some houses in flood risk areas. An illustration for both errors (vertical and horizontal inaccuracy) is shown in Figure 2: on the left the real flood area has been computed, and the image on the upper right shows an error due to vertical inaccuracy (the CompDEM altitude is higher than the RefDEM altitude and an urban area is flooded); the image on the bottom right shows the error due to horizontal inaccuracy (the CompRef is moved to the left and up direction and an urban area is flooded). Another example of horizontal accuracy importance comes from the ability of LiDAR data to produce a highly detailed drainage network [9] and it would be interesting to know the variation in the streams if computed with different DEM sources.
Landslide is another field where horizontal displacement has to be computed. There are several approaches to get such calculations, and one of them is based on series of images taken in different epochs, e.g., [10,11]. Those studies base their computation on minimizing a correlation function between a local window in the reference image and its homologous window in the compared image, according to the method explained by [12].
In our present study we use an approach similar to that of [8], but apply the particle image velocimetry (PIV) algorithm. PIV comes from a distant research field, more commonly used in studies involving aerodynamics [13], transonic flows [14], wind tunnels [15], or aeroacoustic predictions [16]. PIV can be used to determine the velocity field generated by a flow, tracking the particles inside that flow by repeated photography of the same area; an algorithm then computes the displacement of a particle (or a set of particles) from one photograph to the next. Dividing the particle displacement (number of pixels multiplied by the pixel size on the ground) by the time between photos gives us the velocity. Figure 3 shows how two sets of particles are displaced from Photo 1 to Photo 2.
Figure 2. Importance of vertical and horizontal accuracy.
Figure 2. Importance of vertical and horizontal accuracy.
Remotesensing 08 00014 g002
Figure 3. Particle displacement derived from the position captured in two different photographs.
Figure 3. Particle displacement derived from the position captured in two different photographs.
Remotesensing 08 00014 g003
A DEM resembles a digital photograph in two ways:
-
Both of them are made up of a regular grid of cells.
-
Each cell is represented by a numerical value.
For these reasons, the PIV algorithm used on photos can almost directly be applied to DEMs. The only difference is that RGB (Red, Green, Blue) color photos must be transformed to a gray scale in the case of PIV, while DEM does not require transforming the cell values.

2. Materials and Methodology

In order to test the ability of PIV in estimating the horizontal displacement between two DEMs, we adopted the dataset used in [8]. It consists of four samples (corresponding to four different places) pertaining to terrains with different slopes; they are labeled as upland, hill and mountain according to the classification from [17]. Each sample is composed of two DEMs from two different sources: “Instituto Geográfico Nacional” (IGN) [18] and “Instituto Cartográfico de Andalucía” (ICA) [19], Spanish cartographic agencies at the national and regional level, respectively. The data captured correspond to the years 2005 (IGN) and 2004 (ICA). The cell size for the ICA DEMs is 10 m, and for the IGN ones it is 25 m. This means that the IGN DEMs must be resized in order to have the same dimensions in homologous DEMs (we will refer to homologous DEMs to the sample composed for the ICA and the IGN DEMs from the same place). Once we resize an IGN DEM, its cell size is the same as the ICA one, 10 m. There are several interpolating methods that could be used for resizing the IGN DEM, taking into account the cell neighborhood, e.g., the nearest, the bilinear or the bicubic. We decided to use the bicubic one, because the interpolating method influence on the estimated horizontal displacement is small, according to the study shown in Section 3.1.
The names and the areas of each sample DEM are as follows: Seville (63.98 km2), Huétor-Tájar (34.39 km2), Colmenar (70.64 km2) and Quéntar (51.85 km2). Table 1 shows the classification and the slope value for each place.
Table 1. Slopes (%), heights (m) and category for each DEM from the dataset.
Table 1. Slopes (%), heights (m) and category for each DEM from the dataset.
Slope (%)Terrain category UplandHillMountain
Places SevilleHuétor-TájarColmenarQuéntar
Cartographic AgenciesICA (10 m cell size)slope2.668.2417.9323.28
Height50.29566.4367.301330.19
IGN (25 m cell size)slope2.538.417.4122.11
Height50.30566.41367.321330.25
To compare the PIV-estimated horizontal displacement (PIVEHD) with a reference displacement, we compute what we called the “real horizontal displacement” (RHD), which will be considered as the ground truth.
The RHD was computed by digitizing, as a polygonal shape (linear polyline), the same hydrological feature (creek or river) from the two DEM sources (ICA and IGN). Because both polygonals are in the same reference system (ETRS89 Datum) and they do not coincide exactly, the horizontal displacement can be computed as the area comprised between the two homologous polygonals and divided by their mean length. Each hydrological feature was digitized by four different people to enhance RHD accuracy. For the sake of redundancy we compute the RHD for four different hydrological features in each DEM (denoted as hydrological features 1, 2, 3 and 4 in Table 2, distributed as shown in Figure 4). Table 2 shows the RHD between homologous hydrological features and their corresponding length; the last columns give the RHD mean value for each place, taking into account the four hydrological features.
Table 2. Real displacement (RD) between homologous DEMs computed comparing linear hydrological features.
Table 2. Real displacement (RD) between homologous DEMs computed comparing linear hydrological features.
Hydrological Features1234
MeasuresRHD (m)Length (m)RHD (m)Length (m)RHD (m)Length (m)RHD (m)Length (m) M R H D (m)
Quéntar4.456416.924.757574.684.553530.284.012454.974.53
Colmenar5.1410,519.145.365180.194.993537.424.577054.265.01
Huétor-Tájar5.033155.435.671356.235.182151.205.4413,410.005.36
Seville6.845279.106.553782.674.583527.305.0114,224.705.53
Figure 4. Sample DEMs (four different places) and distribution of the digitized linear hydrological features. Surrounding each hydrological feature is a buffer 600 m wide.
Figure 4. Sample DEMs (four different places) and distribution of the digitized linear hydrological features. Surrounding each hydrological feature is a buffer 600 m wide.
Remotesensing 08 00014 g004
It may be argued that applying the PIV to the whole DEM—comparing the PIVEHD of the entire DEM with the RHD from linear hydrological features—is an over-generalization. Hence, we computed the PIVEHD inside a buffer 600 m wide surrounding each hydrological feature. This was compared to the corresponding RHD. The size and shape of the buffers can be seen in Figure 4.
Depending on the type of analysis needed, the PIV algorithm can be set with different parameters. A complete overview of PIV setting variability is found in [20]. We opted to run the algorithm based on the cross-correlation function performed under Fast Fourier Transform (FFT), taking into account that displacement could entail translation (case A in Figure 3) or a mixture of translation and rotation (case B in Figure 3), as explained in [20].
We used PIVlab software [21] to arrive at parameters suitable for the above requirements.
The PIV method consists of defining an interrogation window in the first DEM (D) and displacing it over the second DEM (D′) around the original position of the window in D (Figure 5).
Figure 5. Interrogation window on the first DEM (D) and a displacement over the homologous DEM (D′).
Figure 5. Interrogation window on the first DEM (D) and a displacement over the homologous DEM (D′).
Remotesensing 08 00014 g005
The PIV algorithm attempts to find the best match between images (in our case DEMs D and D′) using the discrete cross-correlation function defined in Equation (1).
R I I ( x ,   y ) = i = K K j = L L H ( i , j ) H ( i + x , j + y )
where H and H′ respectively represent the cell height in D and D′. For each displacement (x, y) selected, we obtain a R I I ( x ,   y ) cross-correlation value. If we select a shift range ( M ) x M , N y N we obtain a correlation plane of size ( 2 M + 1 ) × ( 2 N + 1 ) . The cell value for the correlation plane is the sum of the products of all overlapping cells in D and D′. A graphic illustration is shown in Figure 6. The best match would be the highest cell value in the correlation plane, and therefore the displacement is the corresponding (x, y) shifts that generated such a high cell value.
Figure 6. Cross-correlation plane derived from the cross-correlation function using an interrogation window 4 × 4 in size.
Figure 6. Cross-correlation plane derived from the cross-correlation function using an interrogation window 4 × 4 in size.
Remotesensing 08 00014 g006
The cross-correlation function as expressed in Equation 1 has two weak points:
  • High computational time.
  • It may produce different maximum correlation values for the same degree of matching because the function is not normalized. For example, areas with many (or brighter) particle images will produce much higher correlation values than interrogation windows with fewer (or weaker) particle images.
The solution to the first weak point is to compute the cross-correlation function as a complex conjugate multiplication of their Fourier transforms: R I I = H ^ · H ' ^ * , where H ^ and H ' ^ are the Fourier transforms from H   and H .
The second weak point can be overcome using the cross-correlation coefficient, which normalizes the cross-correlation function and its definition as expressed by Equation (2):
c I I ( x , y ) = C I I ( x , y ) σ H ( x , y ) σ H ' ( x , y )
where
C I I ( x , y ) = i = 0 M j = 0 N [ H ( i , j ) μ H ] [ H ( i + x , j + y ) μ H ( x , y ) ]
σ H ( x , y ) = i = 0 M j = 0 N [ H ( i , j ) μ H ] 2
σ H ' ( x , y ) = i = 0 M j = 0 N [ H ( i , j ) μ H ( x , y ) ] 2
The value μ H is the average of the template and is computed only once while μ H ( x , y ) is the average of H coincident with the template H   at position ( x ,   y ) .
PIVlab [21] software offers the option of applying a multiple pass interrogation window that improves the PIV results in the sense that more particle matching is achieved in the successive passes. Each pass decreases the interrogation window size. We tried three hierarchy passes, where in the next pass the interrogation window was one half of the previous one.

3. Results and Discussion

3.1. Effect of the Interpolation Method for Resampling IGN DEM to ICA DEM on the PIVEHD

In order to study the influence of the resampling method to resize the IGN DEM to the size of the ICA DEM, we compute the differences in height for the two extreme cases, i.e., Quéntar (the steepest DEM) and Seville (the flattest DEM). If we denote a resampled DEM as R x m where x represents the DEM (Quéntar (Q) or Seville (S)) and m represents the interpolation method (nearest (near), bilinear (bil) or bicubic (bic)), we can compute the mean difference in height for a particular DEM and two different methods, although it is most informative if we compute the mean of the absolute value of the height difference. For example considering Quéntar and the bicubic and nearest method, such a mean value will be computed as follows:
μ Q b i c n e a r = m e a n ( a b s ( R Q b i c R Q n e a r ) )
Table 3 shows the possible combination represented by Equation (6).
Table 3. Mean absolute differences between interpolation methods to resample the ICA DEM.
Table 3. Mean absolute differences between interpolation methods to resample the ICA DEM.
μ x b i c b i l (m) μ x b i c n e a r (m) μ x b i l n e a r (m)
Quéntar0.392.782.79
Seville0.060.250.25
A histogram representation referring to Quéntar for the three combinations mentioned in Table 3 is shown in Figure 7.
Figure 7. Histogram from the differences between interpolation methods to resample the ICA DEM, Comparations: (a) Bicubic-Bilinear; (b) Bicubic-Nearest; (c) Bilinear-Nearest.
Figure 7. Histogram from the differences between interpolation methods to resample the ICA DEM, Comparations: (a) Bicubic-Bilinear; (b) Bicubic-Nearest; (c) Bilinear-Nearest.
Remotesensing 08 00014 g007
From Table 3 and Figure 7 we observe than the differences between bicubic and bilinear methods are small and, when the terrain is flat like Seville, the differences between all the methods are also small. For this reason we start to compute PIV displacement for a steep area like Quéntar because if the differences between the methods are small, then in flat areas the differences will be smaller. We compute the PIVEHD using the three different IGN resampled DEMs (nearest, bilinear and bicubic) and compute the mean discrepancy in the absolute value for the module of the vector, i.e., if U m is the matrix containing the horizontal axis component for the displacement vectors, V m is the vertical axis component for the displacement vectors and m is the interpolation method (nearest, bilinear an bicubic), we can compute the effect in computing the PIVEHD using Equation (7), which considers the computation made between the bicubic and nearest interpolation methods:
μ d i f f b i c n e a r = m e a n ( s q r t ( ( U b i c U n e a r ) 2 + ( V b i c V n e a r ) 2 ) )
The parameters used for computing PIVEHD where pass 1:32 × 32; pass 2:16 × 16; pass 3:8 × 8.
The Table 4 shows the values for the differences between methods when computing PIVEHD and the corresponding standard deviation.
Table 4. Mean differences for PIVEHD modules and standard deviations comparing pairs of interpolation methods.
Table 4. Mean differences for PIVEHD modules and standard deviations comparing pairs of interpolation methods.
μ d i f f b i c b i l (m) σ d i f f b i c b i l (m) μ d i f f b i c n e a r (m) σ d i f f b i c n e a r (m) μ d i f f b i l n e a r (m) σ d i f f b i l n e a r (m)
Quéntar0.040.110.160.170.160.18
From Table 4 we observe that the effect on PIVEHD is not significant either computed by the bicubic or bilinear method (4 cm), and just the near method differs with respect to bilinear and bicubic method in the same magnitude (16 cm), which is not a big value considering the cell size for ICA (10 m) and IGN (25 m). So we decide to use the bicubic interpolation method in order to resample the IGN DEM to the ICA DEM size.

3.2. Results Analysis

In order to analyze several interrogation window sizes, we tested three settings classified as High, Medium and Small. The values for such a classification are as follows:
-
High: (pass 1:64 × 64; pass 2:32 × 32; pass 3:16 × 16).
-
Medium: (pass 1:32 × 32; pass 2:16 × 16; pass 3:8 × 8).
-
Small: (pass 1:16 × 16; pass 2:8 × 8; pass 3:4 × 4).
A graphic example of applying the Medium interrogation window to the Colmenar region is shown in Figure 8. The displacement vectors are shown in green.
Table 5 shows the PIVEHD expressed in meters (Medium interrogation window setting) for each place (Quéntar, Colmenar, Huétor-Tajar and Seville), evaluated over each of their hydrological features (1, 2, 3 and 4). Although the PIVlab algorithm outputs are vectors, we compute their modules in order to know the magnitude of the displacement; the comparison shown in Table 5 refers to those modules. Remember, the horizontal displacement has been computed inside the buffer area surrounding each hydrological feature (creeks or river), as shown in Figure 4. Given in Table 5 are the respective RHDs in order to compute the error. The two last columns show the mean absolute (ABS) error for each place and the total mean ABS (error).
Figure 8. Horizontal displacement represented by vectors (green color). Detail from Colmenar DEM.
Figure 8. Horizontal displacement represented by vectors (green color). Detail from Colmenar DEM.
Remotesensing 08 00014 g008
Table 5. PIVEHD from a Medium (pass 1:32 × 32; pass 2:16 × 16; pass 3:8 × 8) interrogation window setting and their discrepancies (ABS (error)) with respect to the RHD.
Table 5. PIVEHD from a Medium (pass 1:32 × 32; pass 2:16 × 16; pass 3:8 × 8) interrogation window setting and their discrepancies (ABS (error)) with respect to the RHD.
1234Mean ABS (Error)TOTAL Mean ABS (Error)
PIVEHD0.920.941.160.95 3.53
QuéntarRHD4.454.754.554.01
Error3.533.813.393.063.45
PIVEHD1.881.671.722.28
ColmenarRHD5.145.364.994.57
Error3.263.693.272.293.13
PIVEHD1.972.451.772.48
Huétor-TajarRHD5.035.675.185.44
Error3.063.223.412.963.16
PIVEHD1.231.181.621.33
SevilleRHD6.846.554.585.01
Error5.615.372.963.684.41
The mean ABS (error) shown in Table 5 seems somewhat high, particularly if compared with the estimations obtained using the optical flow algorithm [8] which had a mean error of 1.13 m. Nevertheless, the optical flow algorithm features a parameter (α) that is similar to a scale factor. For this reason, if we find a suitable scale factor that reduces the errors then PIVEHD can be considered as scalable and a good estimator of displacement. We thus propose to compute the best scale factor using the least squares method, minimizing the function   v T v obtained from A x B = v , where A is the estimation from PIV, B is the RHD, and x is the scale factor. For example, considering Mountain terrain (Quéntar), the corresponding values are as follows: A = [ 0.92   0.94   1.16   0.95 ] T , B = [ 4.45   4.75   4.55   4.01 ] T and v = [ 3.53   3.81   3.39   3.06 ] T
Equation (8) computes x and the standard deviation (s)
x = ( A T A ) 1 A T B
s = v T v d o f
where dof means degrees of freedom.
The scale factors for the High, Medium and Small interrogation windows after applying the least squares methods are shown in Table 6. The values are similar, ranging from 2.937 to 2.388.
Table 6. Scale factor minimizing the PIVED errors.
Table 6. Scale factor minimizing the PIVED errors.
Interrogation WindowHighMediumSmall
Scale Factor2.9912. 9372.388
Standard Deviation1.9671.6771.585
After applying the scale factor to the PIVEHD, the adjusted PIVEHD (APIVEHD) and the corresponding errors are shown (Table 7, Table 8 and Table 9), representing the High, Medium and Small interrogation windows.
Table 7. PIVED after applying the scale factor (APIVED). High interrogation window.
Table 7. PIVED after applying the scale factor (APIVED). High interrogation window.
1234Mean ABS (Error) mTOTAL Mean ABS (Error) mTOTAL Standard Deviation ABS (Error) m
QuéntarAPIVEHD2.352.432.642.57 1.60.79
ABS (error)2.102.321.911.441.94
ColmenarAPIVEHD5.705.424.466.83
ABS (error)−0.56−0.060.53−2.260.85
Huétor-TajarAPIVEHD4.577.944.616.38
ABS (error)0.46−2.270.57−0.941.06
SevilleAPIVEHD2.903.296.253.73
ABS (error)3.943.26−1.671.282.54
Table 8. PIVED after applying the scale factor (APIVED). Medium interrogation window.
Table 8. PIVED after applying the scale factor (APIVED). Medium interrogation window.
1234Mean ABS (Error) mTOTAL Mean ABS (Error) mTOTAL Standard Deviation ABS (Error) m
QuéntarAPIVEHD2.692.763.402.79 1.310.51
ABS (error)1.761.991.151.221.53
ColmenarAPIVEHD5.514.915.046.71
ABS (error)0.370.450.052.140.75
Huétor-TajarAPIVEHD5.777.195.207.30
ABS (error)0.741.520.021.861.04
SevilleAPIVEHD3.603.474.763.90
ABS (error)3.243.080.181.111.9
Table 9. PIVED after applying the scale factor (APIVED). Small interrogation window.
Table 9. PIVED after applying the scale factor (APIVED). Small interrogation window.
1234Mean ABS (Error) mTOTAL Mean ABS (Error) mTOTAL Standard Deviation ABS (Error) m
QuéntarAPIVEHD3.033.063.462.96 1.420.36
ABS(error)1.421.691.091.051.31
ColmenarAPIVEHD4.203.754.065.08
ABS(error)0.941.610.930.511.00
Huétor-TajarAPIVEHD6.857.446.607.84
ABS(error)1.821,771,422.401.85
SevilleAPIVEHD4.614,254,103.92
ABS(error)2.232.300.481.091.53
Comparing Table 5 and Table 8, both corresponding to the Medium interrogation window, we note an improvement (Impv) of 2.22 m through application of the scale factor. This amounts to a substantial improvement (62.9%):
I m p v = e r r o r W i t h o u t S c a F a c t o r e r r o r A f t e r S c a F a c t o r = 3.53 1.31 = 2.22
Given that the cell size of the ICA DEMs is 10 m, and the best error corresponds to the Medium interrogation window (1.31 m), we may conclude that the PIV algorithm is a sound estimator of horizontal displacement between homologous DEMs: the error is around one-seventh of the corresponding cell size. On the other hand, if we consider the IGN cell size, the proportion decreases to one-nineteenth.
In general, no great differences are observed among the errors of the three categories (High, Medium and Small); the dispersion appears to decrease when the interrogation window decreases (see the values for the Total Standard Deviation ABS (error) in Table 7, Table 8 and Table 9). However, if we were to recommend an interrogation window, it would be the Small one, in view of its reduced Total Standard Deviation ABS (error) value. Notwithstanding, the Medium category offers the advantage of a shorter processing time.

3.2.1. Performance in Comparison to Previous Approaches

As mentioned in the abstract and the Introduction section, there are barely any algorithms analyzing horizontal displacement on homologous DEMs as a whole. Just the contour-based [7] and optical flow [8] approaches study the DEM horizontal variable. So we consider it necessary to compare the performance between both algorithms, and such a performance is studied considering the absolute error between the ground truth and the estimation achieved by every approach (contour, optical flow and PIV). An additional comparison concerning computation time has been added in order to assess the global performance. Absolute error and computation time are shown in Table 10. The computer used to test the mentioned algorithms had the following features: Intel i5 processor, 8 Gb RAM, 256 SSD hard drive, Windows 7 operating system.
Table 10. Horizontal displacement absolute errors and computation time obtained using contours, optical flow and PIV approaches.
Table 10. Horizontal displacement absolute errors and computation time obtained using contours, optical flow and PIV approaches.
QuéntarColmenarHuétor-TájarSeville
Contour-basedABS (error)0.520.440.480.52
Computation Time (seconds)548639568677
Optical FlowABS (error)0.390.420.563.18
Computation Time (seconds)122165130151
APIVEHDABS (error)1.530.751.041.9
Computation Time (seconds)143173155184
From Table 10, Mountain and Hill terrains have been similarly estimated by contour and optical flow algorithms, and both get better estimations than the APIVEHD one. On the computation time side, optical flow had the best performance; it is even better than APIVEHD although the differences are not very large (around 14% better than APIVEHD). Although the best horizontal displacement comes from the contour-based algorithm, its computation time is 4.28× larger than optical flow and 3.71× larger than APIVEHD.
Likewise, the error values obtained for Seville, a nearly flat region (2.6% slope), are noteworthy. On such a terrain type, the horizontal displacement errors increase quickly [8]. In our case, the APIVEHD algorithm improves results for nearly flat terrains with respect to the optical flow approach, namely 1.9 m for the APIVEHD algorithm versus 3.18 m for the optical flow algorithm, although the contour algorithm keeps the best performance. We can see from Table 10 that the steeper the terrain the better the optical flow algorithm performs; nevertheless, in near-flat areas the APIVEHD approach achieves the best estimation. Taking into account both absolute error and computation time, we think that our new proposed algorithm APIVEHD gets a better tradeoff when we estimate the horizontal displacement in near-flat areas such as Seville.

3.2.2. Assessing the Influence of DEM Cell Size on the APIVEHD Algorithm Results

After computing the scale factor, the results obtained from the APIVEHD applied to our dataset seem to be suitable for estimating the horizontal displacement (Table 7, Table 8 and Table 9), but there is a variable that could be modified for suitable estimation yield by APIVEHD; this variable is the cell size ratio between the DEMs compared. Our study, resumed in Table 7, Table 8 and Table 9, compared IGN (25 m cell size) and ICA (10 m cell size). In order to evaluate the ability of APIVEHD to get a suitable estimation of horizontal displacement across different cell size ratios, we modified the cell size for ICA to 15 m and 25 m and then we applied the APIVEHD algorithm again, comparing those cell sizes (15 m and 25 m) against the IGN DEM (25 m). According to Section 3.1, the bicubic and bilinear resampling methods are similar and practically do not influence the results yielded by APIVEHD; for this reason, we select the bicubic method in order to compute the derived ICA 15 m and 25 m cell sizes in order to assess the effect of using different cell size ratios between DEMs when computing the horizontal displacement with APIVEHD.
For the sake of readability we just present in Table 11 the horizontal displacement between the following homologous DEMs whose cell sizes are inserted in their names: (1) IGN25m-ICA15m; (2) IGN25m-ICA25m. Note that IGN25m-ICA10m corresponds to Table 8. We selected Table 8 as the reference because the Medium interrogation window has the smaller mean absolute error.
Table 11. APIVEHD obtained by comparison between the homologous pairs DEMs IGN25m-ICA15m and between the homologous IGN25m-ICA25m.
Table 11. APIVEHD obtained by comparison between the homologous pairs DEMs IGN25m-ICA15m and between the homologous IGN25m-ICA25m.
1234Mean ABS (Error) mTOTAL Mean ABS (Error) mTOTAL Standard Deviation ABS (Error) m
IGN25m-ICA15mQuéntarAPIVEHD2.813.633.582.55 1.470.16
ABS (error)1.641.120.971.461.3
ColmenarAPIVEHD4.633.023.255.78
ABS (error)0.512.341.741.211.45
Huétor-TajarAPIVEHD3.777.256.776.78
ABS (error)1.261.581.591.341.44
SevilleAPIVEHD4.934.783.233.31
ABS (error)1.911.771.351.71.68
IGN25m-ICA25mQuéntarAPIVEHD3.213.173.222.64 1.570.22
ABS (error)1.241.581.331.371.38
ColmenarAPIVEHD3.533.292.735.95
ABS (error)1.612.072.261.381.83
Huétor-TajarAPIVEHD3.267.086.983.66
ABS (error)1.771.411.81.781.69
SevilleAPIVEHD7.214.582.863.51
ABS (error)0.371.971.721.51.39
From the comparison of Table 8 and Table 11, we do not observe large variations in the TOTAL mean ABS (error). The error estimating the horizontal displacement when we increase the ICA cell size by 5 m is 1.47, just 0.16 m if compared with original ICA cell size, i.e., the influence of the cell size ratio (IGN25m-ICA10m versus IGN25m-ICA10m) in the horizontal displacement estimation by APIVEHD is just 16 cm. The influence in the case of IGN25m-ICA10m versus IGN25m-ICA25m is 0.26 m. Both effects, 0.16 m and 0.26 m, seem to indicate that the APIVEHD algorithm remains robust to the variation in the cell size ratio.

4. Conclusions

In our study we compute the horizontal displacement between two homologous DEMs coming from different sources. This is important because the horizontal discrepancy in both DEMs puts some derived features in non-corresponding places (e.g., hydrological linear features). Such discrepancies may be estimated by using the PIV algorithm as proposed in this paper. However, a larger number of cases must be computed in order to completely generalize our results. Our findings are a novelty in that an algorithm from the distant disciplines of aerodynamics or water flows is here applied to the DTM field.
The PIV algorithm underestimates the real horizontal displacement, but it does so in a scalar sense, which means that it as a scale factor can be used to arrive at a suitable estimation.
Using PIV and applying the scale factor, the error is around one-seventh of the cell size of the smaller DEMs (the ICA DEMs) and around one-nineteenth in the case of the larger cell size DEMs (the IGN DEMs). We consider these to be acceptable error values.
One remarkable finding is that the PIV approach provides reliable horizontal displacement values for nearly flat terrain types, which are the most difficult ones to compute due to their horizontal displacement.
The recommended interrogation windows are similar to those recommended in the original fields of PIV (aerodynamics, fluid mechanics...), i.e., 32 × 32 [20].
Future works along these lines will involve exploring PIV behavior in the horizontal displacement estimation between DEMs derived from LiDAR (very small cell sizes) and photogrammetry (cell sizes similar to the one used here). It would be interesting to determine whether the three interrogation window categories produce results that agree with the results presented here, or if new dimensions for the interrogation windows are necessary. In the case of LiDAR, because the cell size could be around a few centimeters, it would likely be better to use smaller features than hydrological ones, in order to compute the ground truth for the horizontal displacement, e.g., houses, buildings, fences, small dams, small pieces of roads. Another important point will be the computation time, because if we want to analyze a significant area, the file size will increase considerably, as will the PIV computation. Such a situation could need to parallelize some functions from the PIV software.

Author Contribution

Juan F. Reinoso conceived, designed and performed the experiments; Carlos León and Jesús Mataix did the literature review, preparing the data and analyzed the results.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Gao, J. Resolution and accuracy of terrain representation by grid DEMs at a micro-scale. Int. J. Geogr. Inf. Sci. 1997, 11, 199–212. [Google Scholar] [CrossRef]
  2. Torlegard, K.; Ostman, A.; Lindgren, R. A comparative test of photogrammetrically sampled digital elevation models. Photogrammetria 1986, 41, 1–16. [Google Scholar] [CrossRef]
  3. Aryal, S.K.; Bates, B.C. Effects of catchment discretization on topographic index distributions. J. Hydrol. 2008, 359, 150–163. [Google Scholar] [CrossRef]
  4. Wise, S. Effect of differing DEM creation methods on the results from a hydrological model. Comput. Geosci. 2007, 33, 1351–1365. [Google Scholar] [CrossRef]
  5. McAlister, E.; Domburg, N.; Edwards, T.; Ferrier, B. Hydrological modelling of the river Ythan using Arcinfo Grid. In Proceedings of Sixteenth Annual ESRI User Conference, California, CA, USA, 20–24 May 1996; pp. 253–261.
  6. Reinoso, J.F. A priori horizontal displacement (HD) estimation of hydrological features when versioned DEMs are used. J. Hydrol. 2010, 384, 130–141. [Google Scholar] [CrossRef]
  7. Reinoso, J.F. An algorithm for automatically computing the horizontal shift between homologous contours from DTMs. ISPRS J. Photogramm. Remote Sens. 2011, 66, 272–286. [Google Scholar] [CrossRef]
  8. Reinoso, J.F.; León, C.; Mataix, J. Optical flow algorithm as estimator of horizontal discrepancy between features derived from DEMs: Rivers and creeks as case study. Surv. Rev. 2014, 46, 149–154. [Google Scholar] [CrossRef]
  9. TERC (Tahoe Environmental Research Center). Available online: http://terc.ucdavis.edu/research/storm-water/detention-basin.html (accessed on 24 November 2015).
  10. Bencardino, M. Preparation and processing of optical images (photographs and satellite) for the determination of landslide displacement field. In Proceedings of the Seventh International Conference on Informatics and Urban and Regional Planning INPUT 2012, Cagliary, Italy, 10–12 May 2012; pp. 1202–1214.
  11. Delacourt, C.; Allemand, P.; Casson, B.; Vadon, H. Velocity fields of the “La Clapiere” landslide measured by the correlation of aerial and QuickBirdsatellite images. Geophys. Res. Lett. 2004. [Google Scholar] [CrossRef]
  12. Vadon, H.; Massonnet, D. Earthquake displacement fields mapped by very precise correlation: Complementarity with radar interferometry. In Proceedings of the IEEE 2002 International Geoscience and Remote Sensing Symposium, Toronto, ON, Canada, 24–28 June 2002.
  13. Richard, H.; van der Wall, B.G. Detailed investigation of rotor blade tip vortex in hover condition by 2C and 3C-PIV. In Proceedings of the 32nd European Rotorcraft Forum, Maastricht, The Netherlands, 12–14 September 2006.
  14. Scarano, F.; van Oudheusden, B.W. Planar velocity measurements of a two-dimensional compressible wake flow. Exp. Fluids 2003, 34, 430–441. [Google Scholar] [CrossRef]
  15. Berton, E.; Favier, D.; NsiMba, M.; Maresca, C.; Allain, C. Embedded LDV measurements methods applied to unsteady flow investigation. Exp. Fluids 2001, 30, 102–110. [Google Scholar] [CrossRef]
  16. Schram, C.; Hirschberg, A. Application of vortex sound theory to vortex-pairing noise: Sensitivity to errors in flow data. J. Sound Vib. 2003, 266, 1079–1098. [Google Scholar] [CrossRef]
  17. Li, Z.; Zhu, Q.; Gold, C. Digital TerrainModeling; CRC Press: Boca Raton, FL, USA, 2005. [Google Scholar]
  18. IGN (Instituto Geográfico Nacional de España), 2006. Modelo Digital del Terreno Escala 1:25,000. Ministerio de Fomento. Available online: http://centrodedescargas.cnig.es/CentroDescargas/index.jsp (accessed on 15 July 2014).
  19. De Andalucía, J. Modelo digital del terreno de andalucía. relieve y orografía. ICA [DVD] 2005, 1, 10–27. [Google Scholar]
  20. Raffel, M.; Willert, C.; Wereley, S.; Kompenhans, J. Particle Image Velocimetry. A Practical Guide; Springer: Berlin, Germany, 2007. [Google Scholar]
  21. Thielicke, W.; Stamhuis, E.J. PIVlab-Time-Resolved Digital Particle Image Velocimetry Tool for MATLAB (Version: 1.32). Available online: http://dx.doi.org/10.6084/m9.figshare.1092508 (accessed on 15 April 2014).

Share and Cite

MDPI and ACS Style

Reinoso, J.F.; León, C.; Mataix, J. Estimating Horizontal Displacement between DEMs by Means of Particle Image Velocimetry Techniques. Remote Sens. 2016, 8, 14. https://doi.org/10.3390/rs8010014

AMA Style

Reinoso JF, León C, Mataix J. Estimating Horizontal Displacement between DEMs by Means of Particle Image Velocimetry Techniques. Remote Sensing. 2016; 8(1):14. https://doi.org/10.3390/rs8010014

Chicago/Turabian Style

Reinoso, Juan F., Carlos León, and Jesús Mataix. 2016. "Estimating Horizontal Displacement between DEMs by Means of Particle Image Velocimetry Techniques" Remote Sensing 8, no. 1: 14. https://doi.org/10.3390/rs8010014

APA Style

Reinoso, J. F., León, C., & Mataix, J. (2016). Estimating Horizontal Displacement between DEMs by Means of Particle Image Velocimetry Techniques. Remote Sensing, 8(1), 14. https://doi.org/10.3390/rs8010014

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop