Ability to detect and locate gross errors on DEM matching algorithm

Abstract

Digital elevation model (DEM) matching techniques have been extended to DEM deformation detection by substituting a robust estimator for the least squares estimator, in which terrain changes are treated as gross errors. However, all existing methods only emphasise their deformation detecting ability, and neglect another important aspect: only when the gross error can be detected and located, can this system be useful. This paper employs the gross error judgement matrix as a tool to make an in-depth analysis of this problem. The theoretical analyses and experimental results show that observations in the DEM matching algorithm in real applications have the ability to detect and locate gross errors. Therefore, treating the terrain changes as gross errors is theoretically feasible, allowing real DEM deformations to be detected by employing a surface matching technique.

Keywords:

• DEM matching
• gross error detection
• gross error location
• gross error judgement matrix
• least z-different algorithm
• deformation detection
• robust estimator

1. Introduction

Digital elevation models (DEM) form an important base of the Digital Earth, which was proposed by Al Gore in 1998, and is a powerful tool for studying the hazards on the Earth's solid surface, e.g. Ellen and Mark (1993) used DEMs to map debris flow hazard in Honolulu, Kääb (2000) and Huggel et al. (2003) used DEMs to monitor landslides in the Alps. In recent years, DEM matching techniques (Rosenholm and Torelegård 1988), which can spatially align DEMs without the aid of ground control points, have become a hot topic. Based on this technique, terrain changes can be detected automatically with multi-temporal DEM (Habib et al. 2001, Zhang et al. 2006). This is very important for evaluating the range and degree of motion of the Earth's natural hazards, e.g. landslides and debris flows (NOAA-USGS Debris Flow Task Force 2005).

The widely used algorithm for matching multi-temporal DEMs is the least z-difference (LZD) (Rosenholm and Torelegård 1988), and many extensions of it have been developed for detecting terrain changes. Pilgrim (1996a, b) improved this method by using an M-estimator to replace the least squares (LS) technique used in LZD, naming it the M-LZD method. Li et al. (2001) integrated the LMS-estimator with a random sample scheme and then proposed a new algorithm, called the MS-LZD method. This algorithm can detect nearly 50% of the deformation area. Zhang and Cen (2008) enhanced the performance in both low and high terrain change percentage by using least trimmed square (LTS) and a self-adaptive threshold. To handle terrain changes, all these methods employ one common strategy: taking terrain changes as gross errors, and then replacing the LS technique used in the original LZD with a different robust estimator, e.g. M-estimator, least-of-median-squares (LMS) estimator, or LTS estimator. The deformation detecting ability has been improved greatly. However, according to the surveying errors and reliability theory, without knowledge about how the gross error is detected and located, the results of gross error detection will be unreliable (Li and Yuan 2002). That is to say, without analysing the ability to detect and locate gross errors in the DEM matching algorithm, the detected terrain changes will not be reliable.

There are two possible ways to analyse the detection and locating of gross errors: one is the correlation coefficient (Li and Yuan 2002), and the other is the gross error judgement matrix proposed by Cen et al. (2003). The conclusion of each method is in accordance with the other. As the latter is used before least squares adjustment, it is employed in this paper to study the ability to detect and locate gross errors in the LZD algorithm, a widely used DEM matching algorithm.

The following sections are arranged as follows: Section 2 gives a brief introduction of the coefficient matrix and judgement matrix of DEM matching algorithm; Section 3 presents a theoretic discussion; Section 4 shown two sets of experiments; and some conclusions are given in Section 5.

2. The coefficient matrix and judgement matrix

2.1 Coefficient matrix of LZD

The LZD algorithm uses seven parameters, i.e. three rotation parameters (r _X , r _Y , r _Z ), three translate parameters (t _X , t _Y , t _Z ), and one scale parameter (s), to describe the transformation between DEMs according to the rigid transformation model. Its objective function can be written as

1

where A is the coefficient matrix of the LZD algorithm, V is the vector of residuals, X = [t _X , t _Y , t _Z , r _X , r _Y , r _Z , s] ^T is the vector of the unknown transformation parameters, dz is the vector of constants (Z difference), n is the number of observations. For a real DEM in practice, usually n>>7.

The detailed coefficient matrix can be written as (for more detail, see Rosenholm and Torelegård 1988)

2

where x,y,z are three-dimensional coordinates, and F _X , F _Y are the gradients in the x, y directions, respectively. Every row corresponds to one point.

2.2 The rank of the coefficient matrix

According to Equation (2), the elements of the coefficient matrix A are functions of the plane coordinates, height, and F _X , F _Y . When the DEM is a regular surface, e.g. a plane – the simplest example, the column vectors in the coefficient matrix may be linearly dependent. Then the coefficient matrix A is of deficient rank. In that case, the LZD cannot find the correct matching. Therefore, the ability to detect and locate the gross errors does not need to be discussed. However, a rigorous regular surface does not exist in the real world. Even with a real rigorous regular surface, its corresponding digital surface is still not a rigorous one owning to the inevitable random errors happening during the data collecting and/or digitising processes. Therefore, it is reasonable to suppose that the coefficient matrix A is always full rank, i.e. rank( A ) = 7.

There must exist one maximal linearly independent bundle in A . Suppose that A _t is one of the maximal linearly independent bundles. rank( A _t ) = 7. With the elementary transformation, A _t can be arranged in the seven top rows in the coefficient matrix. Then A can be rewritten as the following block matrix

3

where A _t , A _r are sub-matrices.

A _r is also full rank, rank( A _r ) = 7, for the following two reasons. On the one hand, the column vectors in the coefficient matrix are unlikely to be linearly dependent owing to the complex surface of the Earth, on which the height and gradient change irregularly. On the other hand, the sub-matrix A _r holds n -7 (n>>7) row vectors, only seven row vectors less than A .

2.3 Gross error judgement matrix

For rank( A _t ) = 7, let the inverse matrix of A _t be B

4

where B _i (1≤ i ≤7) is the ith column vector in B, , is the jth element of B _i .

The submatrix A _r can be described as

5

where A _r (8≤i≤n) is the ith row vector in A _r , , is the jth element in A _ri .

According to Cen et al. (2003), the gross error judgement matrix J can be described as

6

where J _k (1≤k≤7) is the kth column vector in J , hereafter called the judgement vector.

Substituting Equation (4) and Equation (5) into Equation (6) leads to

7

3. The ability to detect and locate gross error

3.1 The detecting ability

Proving by contradiction

Suppose that observations in the LZD algorithm do not have the detecting ability, at least one of the following two conditions will be satisfied according to the theory inCen et al. (2003):

1. J is a zero matrix;

2. J _k is a zero vector (1≤k≤7).

First, we disprove condition (2).

Using Equation (7), condition (2) can be represented as

8

For rank ( B ) = 7,

where B _k is the kth column vector of B , is the jth element in B _k .

Considering an arbitrary row a _ri in the coefficient matrix of LZD

9

is made up of the plane coordinates (x,y), height, and F _X , F _Y . The plane coordinates are regularly distributed, but the height and F _X , F _Y are independent, thus

10

So, Equation (8) is not satisfied. Then J _k (1≤k≤7) is not a zero vector. That is to say condition (2) is also not satisfied.

Obviously, J is not a zero matrix, because J _k (1≤k≤7) is not a zero vector. Therefore, condition (1) is also not satisfied.

In summary, neither condition is satisfied. Therefore, the original supposition is false. Observations in the LZD algorithm do have the detecting ability.

▪

3.2 The locating ability

Proving by contradiction.

Suppose that observations in the LZD algorithm do not have the locating ability, at least one of the following three conditions will be satisfied according to the theory in Cen et al. (2003):

1. J is a row vector;

2. ∃ J _k ⊂ J (1≤k≤7) has all of its elements equal to zero except for the i-th (8≤i≤n) element;

3. ∃ J _k , J _m ⊂ J J _k , J _m is linearly dependent (1≤m≠k≤7).

First, according to Equations (6) and (7), J could not be a row vector for n>>7. Therefore, condition (1) cannot be satisfied.

Second, letting the sth (8≤s≤n) element in J _k be , according to Equation (7), we can obtain

11

Then condition (2) can be rewritten as the following two equations:

12

13

The above two equations should be satisfied simultaneously.

Using Equations (11), Equation (12) can be rewritten as

14

The plane coordinate could not equal 0, and also the height and gradient values do not always equal 0, so must not always equal 0, except for

15

However, when Equation (15) is satisfied, (8≤s≤n). This conflicts with Equation (13). Therefore, condition (2) is not satisfied.

Finally, we disprove condition (3).

Because J _k , J _m are linearly dependent, there must exist two real numbers, e, f, (e and f are not equal to zero simultaneously):

16

Substitute Equation (7) into Equation (15)

17

Referring to Section 3.1, we can obtain A _ri ≠0.

Also rank ( B ) =7 ⇒ B _k , B _m is linearly independent ≠0 (1≤m≠k≤7).

Therefore, Equation (17) requires

18

The row vector a _ri is composed of the plane coordinates, height and F _X , F _Y , and the height and F _X , F _Y are independent of each other, so

19

Therefore, Equation (18) is not satisfied. Therefore, condition (3) is also not satisfied.

All three conditions are not satisfied. The original supposition is false. Observations in the LZD algorithm do have locating ability.

▪

4. Experiments and results

4.1 Test I

First of all, we will discuss the above analyses on two simplest data sets. Both of them contain only 3×3 points, which are the least size of the surface LZD algorithm can deal with. The height value of each point is shown in Figure 1, and the gradient of each point in X, Y directory is listed in Table 1. For convenience, all points in both DEMs are numbered from 1# to 9#.

Figure 1.  Data set.

Table 1. Gradient of DEM.

DEM	Gradient in X directory			Gradient in Y directory
Natural DEM	−1.01	−0.25	−0.11	0.61	0.64	0.85
	−0.98	−0.31	−0.26	−0.25	−0.18	−1.05
	−0.11	−1.04	−1.29	−0.31	−0.41	−1.58
Plane DEM	−0.21	−0.21	−0.21	0.20	0.20	0.20
	−0.21	−0.21	−0.21	0.20	0.20	0.20
	−0.21	−0.21	−0.21	0.20	0.20	0.20

The left-hand side of Figure 1 is extracted from real DEM, called Natural DEM and the right one is a simulated strict inclined plane, called Plane DEM. They are used as the reference DEM, and the ready matching DEM is derived by applying the given transformation parameters, and then Gauss distribution random errors are added to the reference DEM.

According to Equation (6), their judgement matrix J can be derived, and listed in Tables 2 and 3.

Table 2. The judgement matrix of natural DEM.

	Essential observations
Redundant observations	1#	2#	3#	4#	5#	6#	7#
8#	0.3	−2.7	1.3	0.5	3.5	−0.4	−1.5
9#	1.5	−6.6	2.8	−0.5	5.5	0.5	−2.2

Table 3. The judgement matrix of plane DEM.

	Essential observations
Redundant observations	1#	2#	3#	4#	5#	6#	7#
8#	0.0	1.0	0.0	0.0	−2.0	0.0	0.0
9#	−1.0	2.0	0.0	2.0	−4.0	0.0	0.0

There are no zero elements in the judgement matrix of natural DEM, and some column vectors are linear dependent. There are three zero-column vectors in the judgement matrix of plane DEM, and the gross error in their corresponding observations could not be detected.

To discuss this in depth, the LZD's coefficient matrix A on plane DEM is given in Equation (20)

20

Obviously, the 1^st, 2^nd, 3^rd column vectors are linear dependent on each other, and the 4^th, 5^th and 6^th column vectors are also linear dependent (A ₄ – 0.20A ₅ + 0.21 A ₆ = 0). Therefore, the coefficient matrix is of rank defect, i.e. the number column rank is 4 and the number of rank defect is 3. This is because the gradient in the X, Y directory of plane DEM is constant. According to the above section, if there is no zero-column vector in judgement matrix, the coefficient matrix should be full rank.

From the coefficient matrix A (Equation (16)), the rank defend only exists in a rigorous surface; the plane DEM used in the above tests is a typical example. The gradient in X, Y directory is varied regularly with its regular height.

To discuss the effect of the random errors on the ability to detect and locate gross errors on the LZD algorithm, the plane DEM used in the above test are taken as the data set. First, the reference DEM is derived by adding zero-mean Gauss-distribution random errors to the plane DEM. The random errors are subject to N ~ (0.0.15²m²). Then the gradient of the plane DEM with errors are listed in Table 4.

Table 4. The gradient of plane DEM with random errors.

DEM	The gradient in X direction			The gradient in Y direction
plane DEM with random errors	−0.22	−0.23	−0.24	0.19	0.21	0.22
	−0.21	−0.21	−0.22	0.19	0.20	0.22
	−0.20	−0.20	−0.20	0.20	0.20	0.21

Compared with the gradient listed in Table 1 and Table 4, the gradient is not constant. Although the relative standard error (RSE) of added random errors is less than that of related criterion (Li and Zhu, 2003), the linear dependent between column vectors in coefficient matrix A is destroyed, the coefficient matrix is of full rank. The new judgement matrix is listed in Table 5.

Table 5. The new judgement matrix.

	Essential observations
Redundant observations	1#	2#	3#	4#	5#	6#	7#
8#	33.0	41.3	−67.3	−73.7	26.0	34.7	7.0
9#	45.9	58.7	−95.6	−103.5	37.4	49.3	8.8

As shown in Table 5, there are no zero elements in the new judgement matrix, and arbitrary two column vectors are linear independent. The new judgement matrix illustrates that a gross error in an arbitrary point could be detected, and can locate one gross error.

4.2 Test II

To give a full and in-depth analysis, four DEMs (Figure 2), containing three typical terrain categories, i.e. slope, valley, and ridge, are adopted to perform more experiments. Their parameters are listed in Table 6. All of them are extracted from a real DEM derived from the aerial photographs in PUWAIGOU (for more detail, see Zhang et al. 2006) in 1957, and they acted as the reference DEM. The ready matching DEM is also sourced from aerial photographs at the same site in 1987. During the past 30 years, debris-flow has erupted many times, and terrain surface has changed greatly.

Figure 2.  Experimental DEM.

Table 6. The parameters of DEM.

DEM	Maximum height /m	Minimum height /m	Average gradient /degree	class
A	2750	2450	26.2	ridge
B	3070	2270	27.3	ridge
C	2450	2230	23.6	valley
D	2720	2360	18.8	slope

The adopted four DEMs have uniform size, 117 by 100, with uniform interval equal to 10 m. According to the character of LZD, the marginal points are dropped during the matching process. Then the actual number of data points used is 115×98=11270.

In this group test, we analyse the ability to detect and locate gross errors on a DEM matching algorithm using its coefficient matrix, and then detect the terrain deformation associated by the LTS estimator to verify whether the deformation can be detected or not, and the detected deformation is corrected or not. To judge the detected deformation, about 2000 check-points are randomly selected, similar to that in Zhang et al. (2006). All of the experimental results show the error of detected deformation is about equal to the random error of DEM used. The analysis of ability to detect and locate gross error is given below in detail.

In each trial, seven evenly distributed arbitrary points are selected as essential observations, and then the judgement matrix is calculated. To each DEM, about 30 trials are performed. One judgement matrix on DEM A is given in Table 7.

Table 7. The judgement matrix.

	Essential observations (Point #)
Redundant observations (Point #)	2040#	2070#	2100#	4500#	8100#	8130#	8160#
1#	13.69	−10.73	−3.14	−2.64	3.23	−1.77	−0.28
2#	11.20	−8.56	−2.57	1.56	2.58	−1.34	−0.30
3#	9.96	−7.51	−2.38	0.75	2.28	−0.98	−0.38
…	…	…	…	…	…	…	…
11268#	−17.62	15.01	4.31	3.56	−3.32	1.78	0.84
11269#	−18.18	15.47	4.44	4.21	−3.42	1.79	0.90
11270#	−16.96	14.38	4.16	2.63	−3.05	1.68	0.79

All the trials gave similar results. The judgement matrix is not zero matrix or row vector, and does not have zero column vector. Also arbitrary two column vectors are linear independent. This result shows that the gross error in arbitrary point is detectable and locatable, and is consistent with the conclusion drawn in the above section.

5. Conclusion

In this paper, the gross error judgement matrix is employed to give an in-depth analysis about the ability to detect and locate gross error in the LZD algorithm. Both the theoretic prove and the results on simulated and real data sets are illustrated LZD algorithm has the ability to detect and locate gross error on real terrain DEM.

The coefficient matrix of LZD on real DEM is sure of full rank, and no zero elements in judgement matrix. Considering the number of gross errors is usually no more than 50%, each gross error is both detected and located in real application. That is to say, if the gross errors can be detected at prearranged confidence levels by the robust estimator, it can be identified correctly.

DEM matching algorithm using least squares could be extended to detect terrain deformation by robust estimator between multi-temporal DEMs, and the detected deformation is reliable.

Notes on contributors

Tonggang Zhang obtained dual BS degrees in Photogrammetry and Remote Sensing and Computer Science and Technology (2000), and then the PhD degree in Geodesy and Surveying Engineering (2006) from SWJTU. He has been a lecturer of the Department of Surveying Engineering, Southwest Jiaotong University, China from 2006. His research interest includes robust DEM co-registration, deformation detection, LIDAR, and pattern recognition.

Minyi Cen received BS degree in Aerial Photogrammetry, MS degree in Engineering Surveying and PhD degree in Road and Railway Engineering from the Southwest Jiaotong University, China, in 1982, 1990 and 1999, respectively. He is a professor of the Department of Surveying Engineering, Southwest Jiaotong University. He worked in the Hong Kong Polytechnic University as a research associate from 2000 to 2001. His research interest includes Error theory, robust image or DEM matching, GPS application, Accurate engineering surveying.

Zizhen Ren obtained M.S. degree in Cartography and GIS Engineering from Southwest Jiaotong University in 2005. She is currently a PhD candidate. Her research work mainly focusew on the data processing of LiDAR.

Ronghao Yang obtained BS degree in Surveying and Mapping Engineering from Chengdu University of Technology, Chengdu, China, in 2002 and the M.S. degree in Computer Application Technology from Southwest Jiaotong University in 2005. He is currently a Ph.D. candidate of Southwest Jiaotong University and a lecturer of the Department of Surveying and Mapping, Chengdu University of Technology, China. His current research interests are in digital photogrammetry and pattern recognition, artificial intelligence and application of engineering surveying technology in geological disaster monitoring and early warning.

Yicong Feng obtained MS degree in Cartography and GIS Engineering from Southwest Jiaotong University, China, in 2003. He has been, since 2005, an Engineer in Land and Resource Department of Sichuan Province, China. He has been, since 2007, a PhD candidate in the Department of Surveying Engineering, Southwest Jiaotong University, China. His research work mainly focuses on theory and application of geographic information system (GIS), DEM matching and its applications, especially on electronic government (E-Government) based on GIS.

Jun Zhu obtained PhD degree in Cartography and the GIS from the Graduate of the Chinese Academic of Sciences in 2006, and M.S. degree from Southwest Jiaotong University, China, in 2003. In August 2007 to July 2008, he was a Postdoctoral Research Fellow in the Institute of Space and Earth Information, Chinese University of Hong Kong. He has been, since 2006, a lecturer in the Department of Surveying Engineering, Southwest Jiaotong University, China. His research work mainly focuses on theory and technology of three-dimensional geographic information system (3DGIS), especially on collaborative virtual geographic environment (CVGE). Dr. Zhu organises and participates in more than 10 scientific and technological projects, and has published 26 research articles by first author.

Acknowledgement

This research is supported by the National High Technology Plan (863) of the People's Republic of China, Project No. 2009AA12Z207.