Show Summary Details

Page of

Printed from Oxford Research Encyclopedias, Natural Hazard Science. Under the terms of the licence agreement, an individual user may print out a single article for personal use (for details see Privacy Policy and Legal Notice).

date: 10 April 2021

Persistent Homology for Land Cover Change Detectionfree

  • Djamel BouchaffraDjamel BouchaffraCentre for Development of Advanced Technologies
  •  and Faycal YkhlefFaycal YkhlefCenter for development of Advanced Technologies, Algeria


The need for environmental protection, monitoring, and security is increasing, and land cover change detection (LCCD) can aid in the valuation of burned areas, the study of shifting cultivation, the monitoring of pollution, the assessment of deforestation, and the analysis of desertification, urban growth, and climate change. Because of the imminent need and the availability of data repositories, numerous mathematical models have been devised for change detection. Given a sample of remotely sensed images from the same region acquired at different dates, the models investigate if a region has undergone change. Even if there is no substantial advantage to using pixel-based classification over object-based classification, a pixel-based change detection approach is often adopted. A pixel can encompass a large region, and it is imperative to determine whether this pixel (input) has changed or not. A changed image is compared to the available ground truth image for pixel-based performance evaluation. Some existing change detection systems do not take into account reversible changes due to seasonal weather effects. In other words, when snow falls in a region, the land cover is not considered as a change because it is seasonal (reversible).

Some approaches exploit time series of Landsat images, which are based on the Normalized Difference Vegetation Index technique. Others evaluate built-up expansion to assess urban morphology changes using an unsupervised approach that relies on labels clustering. Change detection methods have also been applied to the field of disaster management using object-oriented image classification. Some methodologies are based on spectral mixture analysis. Other techniques invoke a similarity measure based on the evolution of the local statistics of the image between two dates for vegetation LCCD. Probabilistic approaches based on maximum entropy have been applied to vegetation and forest areas, such as Hustai National Park in Mongolia. Researchers in this field have proposed an LCCD scheme based on a feed-forward neural network using backpropagation for training. This paper invokes the new concept of homology theory, a subfield of algebraic topology. Homology theory is incorporated within a Structural Hidden Markov Model.

Previous Work and Motivation

This article addresses the question of whether irreversible change of land cover can be detected through construction of structural information using multidate set of images (Figure 1). Several studies have been conducted to map and identify landscape changes. Volpi, Tuia, Bovolo, Kanevski, and Bruzzone (2013) addressed the problem of high-resolution ground cover change detection using support vector machine classification with radial basis function. Lu, Mausel, Batistella, and Moran (2005) developed a change detector based on a combination of per-pixel impervious surface mapping with filtering and unsupervised classification. Du, Xia, and Li (2015) introduced the subpixel aspect in change detection for urban land cover analysis using a back-propagation neural network. Gong, Su, Jia, and Chen (2014) covered change detection systems using synthetic aperture radar images. The latter approach is based on fuzzy c-means clustering using a novel Markov random field energy function. Lv, Zhong, Zhao, Jiao, and Zhang (2016) proposed an accurate detection based on a conditional random field model using high spatial resolution (HSR) imagery. Zhong, Liu, Zhao, and Zhang (2015) proposed pulse-coupled neural networks and the normalized moment of inertia feature dedicated to HSR remote sensing imagery. More recently, Gong, Zhao, and Liu (2016) proposed an LCCD system based on deep neural networks. It is worth noting that Deep Learning (DL) has become the preferred algorithm in data analysis and prediction. A key benefit of DL is its ability to model data through different levels of perception (stacked layers of transformations; Goodfellow, Bengio, & Courteville, 2016). However, DL requires an important number of samples in the training phase to build an efficient classification model. In the context of change detection, data remain scarce and therefore a DL-based model may not outperform a shallow learner.

Figure 1. Land cover change detection using multidate images.

Source: United States Geological Survey (USGS)

This article introduces a novel formalism that merges probability theory (applied mathematics) with topology (pure mathematics). This approach incorporates homology theory within a structural hidden Markov model framework for land cover change detection in remote sensing. This task highlights both the effect of desertification caused by human activities as well as climate change in land cover variation. The proposed probabilistic approach, the topological hidden Markov model, predicts ordered sequences of pixels by providing clues about the shapes they form. Shape disclosure is achieved through a Vietoris-Rips (VR) simplicial complex constructor. The latent local structures are unfolded gradually throughout the observation sequence via these simplicial complexes. These latter polytopes are characterized topologically (via a persistent homology) by their Betti numbers. Experimental results demonstrate the efficiency of the proposed methodology compared to other state-of-the-art land cover change detectors.

The proposed approach to change detection fits within the realm of machine learning supplemented by clues from homology theory (Prasolov, 2007). It is different from most of the traditional available techniques in that radiometric structural changes are captured using topological qualitative features. The focus is on the search for these shapes (called local structures) formed by the data via the exploration of topological features predominant on airborne maps (Figure 2). Topological attributes are extracted via (a) VR simplicial complexes construction; (b) the concept of homology, which measures their Betti numbers (De-Alarcon, Pascual-Montana, Gupta, & Carazo, 2002; Edelsbrunner & Harer, 2010; Hajij, Wang, & Scheidegger, 2018); and (c) the persistent homology captured by cross-validation. Wavelet decomposition is performed to generate a cloud of points assigned to a pixel.

A novel mathematical framework is proposed that locally incorporates topological clues (or signatures) within a two-level HMM framework (Bengio & Frasconi, 1996; Bouchaffra, 2010, 2012; Bouchaffra, Govindaraju, & Srihari, 1999; Fine, Singer, & Tishby, 1998; Rabiner & Juang, 1993; Sanches, 2000). This mathematical approach handles the data in order to locally disclose their hidden persistent structures. The main contribution of the proposed work is fourfold; it consists of (a) generating a cloud of points through a discrete Wavelet decomposition, (b) assigning VR simplicial complexes to this cloud of points (observations), (c) constructing local structures via the use of persistent homology, and (d) building a structured HMM that integrates local structures and predicts whether a visible observation sequence has changed or not (Figure 3).

Figure 2. Map of land cover shape areas that depicts the presence of geometric and topological features.

Topological Background

Some basic topological definitions deemed necessary to comprehend the subsequent work are provided in this section. The reader is advised to refer to Edelsbrunner (2014) for more comprehensive details.

Simplicial Complexes

Definition 1: A set of k+1 points, x0,x1,,xk, is affinely independent if the k vectors, x1x0,x2x0,,xkx0, are linearly independent.

Definition 2: A k-simplex σ is the convex hull of k+1 affinely independent points. The value of k is the dimension of the k-simplex σ and the points x0 to xk are called vertices.

Simplices of dimension 0,1,2,3 correspond to vertices, edges, triangles, tetrahedral, respectively.

Definition 3: A face of σ is a simplex spanned by a subset of the vertices of σ.

Definition 4: A simplicial complex is a finite collection of simplices, K, such that:


every simplex σK, every face of σK;


for every two simplices σ, ϕK, the intersection, σϕ, is either empty or a face of both simplices.

Figure 3. The entire set of modules involved in the proposed approach.

The dimension of K is the largest dimension of any simplex that is an element of K. The Euler characteristic of K is the alternating sum of simplex numbers:


where k denotes the dimension of K and si corresponds to the number of i-simplices, with 0ik.

Abstract Simplicial Complex

Assume a finite collection of elements called vertices.

Definition 5: An abstract simplicial complex is a system of subcollection, A, such that uA and vuvA. The sets u are called abstract simplices.

The dimension of an abstract simplex is (card-1), where “card” represents the simplex cardinality. However, the dimension of A corresponds to the maximum dimension of any of its abstract simplices. The simplicial complex A is drawn in a Euclidean space via a mapping of (a) each vertex to a point in the d-dimensional Euclidean space d and (b) each abstract k-simplex to the convex hull of the k+1 corresponding points. If this drawing obeys the two conditions of definition 4, then it is referred to as a geometric realization of A. For example, a one-dimensional abstract simplicial complex is a graph, and a geometric realization is its straight-line embedding.

Vietoris-Rips Simplicial Complex

The geometrical constructors known as VR simplicial complex are defined as follows.

Definition 6: Letting R0 be a real number, and S a finite set of points in 2, the VR complex of S and R, denoted as Vietoris-Rips(R), consists of all abstract simplices in 2S whose vertices are at most a distance 2R from one another.

In other words, any two vertices at a distance at most 2R from each other are connected by an edge, and a triangle (or higher-dimensional simplex) is addended to the complex if all its edges are in the complex (Figure 4).

Figure 4. The Vietoris-Rips complex associated to six points equidistant on the unit circle.

Persistent Homology and Betti Numbers Computation

The topological characterization of simplicial complexes is performed through the concept of persistent homology whereby Betti numbers are computed. The theory of persistence is inspired by Morse theory and associated with spectral sequences. It brings a quantitative constituent to the Betti numbers that is steady under disturbances. It aims to disclose the different scales under which a cloud of points can be perceived in one and only one coherent formalism.

Let C be a simplicial complex such as (a) isolated point, (b) line segment, (c) triangle, and (d) tetrahedron obtained using VR constructor. Larger simplicial complexes are made up of these basic ones using different values of d (Figure 5). Therefore, several simplicial complexes are exhibited in order to approximate the true shape of the cloud of points (data set). The ultimate goal of persistent homology consists of determining the one that is the most representative of the true shape formed by the cloud of points. It corresponds to the one that produces the qualitative noiseless features.

Figure 5. 2D simplicial complexes constructed using five points with different values of d. Their building blocks are: (a) points, (b) line segments, and (c) triangles. Betti numbers for each simplicial complex are depicted. This drawing has been generated via MATLAB JAVAPLEX (Homology) ToolBox.

Birth and Death

The concept of birth and death is addressed in this section. The starting point is a sequence of homomorphisms connecting the homology groups of the complexes in the filtration. For ij, Ki represents a subcomplex of Kj, which can be viewed as an injective map fi,j:Ki_Kj. It maps a simplex in Ki to the same simplex in Kj. This is why it is called the inclusion map. It carries over to an inclusion map on the p-cycles, fi,jp:Zp(Ki)_Zp(Kj). This in turn creates a map on homology: ϕi,jp:Hp(Ki)Hp(Kj), which does not often represent an inclusion map. More precisely, if δ is a class in Hp(Ki), and zZp(Ki) is a representative cycle, let ϕi,jp(δ) be the class in Hp(Kj) that contains fi,jp(z). It accepts as input a p-cycle in Ki and sends it to Kj For example, if δ encloses a hole in Ki that fills up at the time Kj is reached, then ϕi,jp maps δ to 0Hp(Kj). The image of ϕi,j p is called a persistent homology group. It includes all p-dimensional homology classes that already have representatives in Ki. A class δHp(Ki) is born at Ki if δ is not present in the image of ϕi1,ip, and a class born at Ki dies entering Kj+1 if ϕi,jp(δ) is absent in the image of ϕi1,jp but ϕi,j+1p(δ) is present in the image of ϕi1,j+1p.

Betti Numbers

A simplicial complex exhibits a group structure through the addition of k-simplices. The resulting free group is called the chain group, Ck. The notion of Betti numbers emanates from the quotient group Hk=Zk/Bk, in which Zk denotes the group structure of all k-chains that have empty boundary; and Bk, subgroup of Ck, is also a subgroup of Zk, called the boundary group. The number of distinct equivalent classes of Hk corresponds to the k-th Betti number denoted by βk. Indeed, the k-th Betti number tallies the number of k-dimensional holes in the simplicial complex. In the case where k=0, β0 counts the number of path-connected components of the simplicial complex. In the case of subsets of the Euclidean space 3, β1 counts the number of independent tunnels (holes) and β2 counts the number of cavities, or enclosed voids. For example, in the case of the solid torus, the Betti numbers are β0=1, β1=1, and β2=0. Therefore, one can state that a torus is topologically represented by the series 1,1,0. However, the surface of a torus has β0=1 (one connected component), β1=2 (the one in the center and the one in the middle), and β2=1 (the inside of the donut). The idea of describing local structures via Betti numbers is much more efficient than the one using Euler characteristic (a connectivity measure that is a zero-dimensional Minkowski functional).


Betti numbers are computed for each simplicial complex (by varying d). From each stage to the next, pairing up the births and deaths, a set of intervals or bars called the barcode of the filtration is obtained. Each bar represents a class in one of the homology groups and thus has a definite dimension. Figure 6 depicts an example of a cloud of five points representing a “house with 1 hole.” The persistent homology (barcodes) discloses for 2d2.8, one connected component (β0=1), one hole (β1=1), and no “voids” or “cavities” (β2=0), yielding a series of Betti numbers 1,1,0.

Figure 6. Chain of simplicial complexes with their barcodes. The persistent homology discloses one connected component (β0=1), 1 hole (β1=1).

Local Structures and Persistent Homology

Because the mission consists of identifying change in the local structures assigned to sequences of observations when viewed at two different points in time, it is therefore crucial to characterize a local structure mathematically. The question one investigates is twofold: (a) what would a local latent structure be and (b) when would this local structure be obtained when there are many possible d-values (multiple resolutions) associated with a VR simplicial complex? Since a latent local structure corresponds to a conformation made of a group of symbols depicting a visible observation sequence, it is reasonable to invoke a geometric constructor (such as VR simplicial complex) to infer the true nature of this shape. All basic symbols of an observation sequence are therefore viewed as a cloud of n-dimensional points in a Euclidean space. The simplicial complexes are built from this cloud and then qualitative topological information (Betti numbers) is successively extracted from them. However, the optimal d* value (that minimizes the error rate) among all possible d-values has been computed using a cross-validation scheme. This choice enables extraction of the most informative qualitative features, which are merely the persistent ones across multiple resolution levels d.

Proposition: A feature vector that contains persistent topologic signatures represents an element of a cluster: This cluster regroups similar simplicial complexes (Betti numbers).

For each visible sequence O=O1O2OT, T simplicial complexes are thus constructed; each one is associated to a visible observation sequence Oi(i=1,,T). For a fixed d value, VR simplicial complexes are finally captured and represented by feature vectors (made of topological attributes) that are further partitioned using a Dirichlet mixtures distribution clustering algorithm (Gorur & Rasmussen, 2010).

Figure 7. Each point is a feature vector of Betti numbers (β0,β1,β2) assigned to a simplicial complex. Each cluster Si in this grouping scheme is a local structure. Three local structures—S1, S2, and S3—are depicted.

Definition 7: A local structure assigned to a simplicial complex (or to an observation sequence Oi) is a cluster Si of feature vectors produced by the optimal DMD clustering.

Once the optimal clustering is computed, each cluster Si (grouping similar shapes) is defined as a local latent structure. Each shape is assigned to a constituent Oi of the entire observation sequence O (Figure 7).

Topology within a Hidden Markov Model Framework

To better position this contribution, the traditional concept of HMM is briefly introduced as a benchmark to the proposed methodology.

Traditional HMM

The HMM approach views an observation sequence O=o1o2o3oT as a series of symbols produced by a sequence of hidden states with a certain emission probability distribution. Each symbol is emitted from one single hidden state using this probability distribution, and the hidden state distribution is a Markov chain. The following definition exhibits the statistical parameters of a tradition HMM.

Definition 8: An HMM is a triplet characterized by λ=[π,A,B], where

π=[πi] is the initial state probability vector,


A=[aij] is the hidden state transition probability matrix from state qi to state qj,1<i,j<N,


N is the number of hidden states in the model, and

B=[bj(t)] is the state conditional probability matrix of the emitted symbol ot given a hidden state qj, with (tbj(t)=1)j, T is the number of observations oi, each one of length ri.

As an example, Figure 8 depicts a simple HMM graph.

Figure 8. A simple HMM Graph with two hidden states and three symbols.

Three problems are assigned to an HMM: (a) probability evaluation, (b) statistical decoding, and (c) parameter estimation (or training). In summary, the evaluation problem in the HMM consists of evaluating the probability for the model λ=[π,A,B] that most likely produces the sequence O. This probability can be expressed as


Using state conditional independence of the observation sequence O, one obtains


The evaluation problem is based on the observation state conditional independence; however, there are numerous applications in which the conditional independence condition does not hold.

Topological HMM: Mathematical Framework

Because the proposed approach extracts qualitative features (using homology theory) from local shapes described by simplicial complexes, it is called a topological hidden Markov model (THMM). To provide a complete introduction of the THMM formalism, this section discusses the THMM context, background on the VR simplicial complexes constructor, and search and characterization of local structures.

THMM Context

In order to capture shapes, the THMM approach views an observation sequence O=O1O2OT of length T as a series of constituents Oi made of strings of symbols oij interrelated in some way. Each observation sequence O is not only one sequence in which all constituents are conditionally independent but also a sequence that is divided into a series of T strings Oi(Oi=oi1oi2oiri), (1iT), (ri is the length of Oi). Each Oi is produced from a set of local structures Si. The observed symbols oij are emitted from hidden states qij.

Sequence Evaluation

Let O=(O1O2OT) be the entire sequence of observations. Each observation Oi is by itself a sequence of symbols oij, (i.e.,Oi=oi1oi2oiri) and is assigned to a structure Si. In other words, the entire visible observation sequence O reveals a sequence of local structures S=S1S2ST. The THMM takes into account the observation sequences with their latent structures. Therefore, the classification problem can be cast into the determination of the class ωi (i=1,2,,c, where c is the total number of classes) that maximizes the probability of the entire observation O given a class model ω. This statement can be written as


However, because (a) of the conditional independence of the Oi with respect to the event {S,ω} and (b) each Oi depends only on its corresponding local structure Si, Prob(O|S,ω) and Prob(S|ω) are computed as


Likewise, it is reasonable to assume that the local structure distribution is a Markov chain of order 1, and therefore one can write


Finally, by regrouping the product terms in equations (2) and (3), one can express Prob(O,S|ω) as follows:


It may be helpful to focus on the different terms separately. It is worth underscoring that Prob(Oi|Si) is difficult to estimate; however, the converse Prob(Si|Oi) conveys more meaning because it invokes properties of the VR constructor. Therefore, using Bayes’s rule, one can write

Prob(Oi|Si)=Prob(Si|Oi)×Prob(Oi)Prob(Si),in whichProb(Si)is a constant (normalizing term).(8)

Therefore, equation (6) can be transformed into


One can decompose Prob(Oi) as


where Oi=oi1,oi2,,oiri, and Qi=qi1,qi2,,qiri. Therefore, equation (10) can be decomposed into


It is assumed that an observation oij depends only on the hidden state qij from which this observation has been emitted. Likewise, a hidden state qij depends only on the previous hidden state qi(j1) (first-order Markov chain). Using these assumptions, equation (11) is transformed into the following:


This latter expression requires an estimation of Prob(o1|qi) and Prob(qi|qi1). Finally, one derives the classification criterion, which is


Equation (12) allows classifying test observation sequences after an optimal training is conducted. Figure 9 depicts the graph of a THMM.

Definition 8: A topological HMM is a quintuplet γ=[π,A,B,C,D] such that

vector π=[πi]=[prob(q1),,Prob(qn)]T; matrixA=[aij]=Prob(qi|qj); matrixB=[bij]=Prob(oi|qj); matrixC=[cij]=Prob(Sj|Si); matrixD=[dij]=Prob(Si|Oj) with the natural constraints as in the traditional HMM framework.

Figure 9. Graph associated with the topological-based hidden Markov model.

Parameter Estimation

The parameters π, A, and B are computed in the same way as in traditional HMM using a Baum-Welch parameter estimation (or expectation maximization), however the matrices C=[cij] and D=[dij] that unfold the local structures are estimated as follows:

d˜ij=number of times knearest neighbors associated to Ojcontained in Si|cluster Si|,

c˜ij=number of occurences of the pair<Si,Sj>inalltrainingsequences<Si,+>,

where |cluster Si| designates the cardinality of cluster Si and <Si,+> represents the number of times the cluster Si occurs in the training set.

THMM Change Detection

This section explains the phases needed to build the THMM for land cover change in remote sensing.

Remote Sensing Feature Generation

Before any classification task such as change detection can take place, one has to extract discriminative and less costly features. It is well known that the classifier power generalization depends in part on the dimension of the feature space. Feature generation (aka extraction) in change detection represents a challenging and tedious task due to diverse regional characteristics such as (a) area scattering: urban regions usually account for small areas and are scattered in different locations; (b) large regions: vegetation and bare soil region represent large and homogeneous surfaces; and (c) thematic similarity: spectral features between some of the thematic classes such as “less dense urban” and “bare soil” are very much alike especially when the spatial resolution is limited. In traditional approaches, feature extraction in change detection is usually performed using software such as ENVI and eCongnition. These programs invoke spatial features (e.g., texture, area, orientation, size) as well as spectral attributes (e.g., histogram depicting occurrence frequencies of gray values, pixel radiometric values in different bands). Because this work is driven by the extraction of local structures that account for image spatial changes between two points in time, it is deemed sufficient to rely on radiometric values at different bandwidths to represent a pixel in an image.

Building VR Simplicial Complexes and Local Structures from a Pixel

The amount of variation of pixels’ radiometric values is an indicator of the type of change in the image. For example, in the case of a disaster, the image exhibits a drastic radiometric change, whereas in the case of urban growth, this radiometric variation is rather smooth. Each image pixel contributes to the formation of a simplicial complex. An observation sequence that represents a series of radiometric vectors in three bands is associated with a pixel. This sequence represents a cloud of points from which a simplicial complex is built. A two-dimensional discrete wavelet transform (DWT) is applied to the entire image of dimension (r,c) for each band out of three. Therefore, this decomposition produces an average image and a sequence of detailed images in each scale. The inverse function of the wavelet transform is applied to these images in order to restore the original dimension (r,c). The goal of this decomposition is to associate each pixel of the original image to its corresponding pixel in a one-to-one mapping.

Let D=[Aj,{Djk},j=1,,J and k=1,,K] be the decomposition in which Aj refers to the approximation image at the j-th scale and Djk is the detail image at scale j and orientation k. For the work described herein, a wavelet decomposition with four levels (scaleJ=4) and three orientations (K=3) for each level, was undertaken. For each scale, the i-th pixel in the original image is associated with a vector containing its three values in the detail images. This procedure is repeated for each of the four scales and then for each band. Hence, DWT decomposition produces a set of 12 points in the details vector space. Figure 10 shows an example of a simplicial complex (built from the 12 points generated from a single pixel using DWT) in the 3D frequency bands. Betti numbers’ qualitative features are extracted from this simplicial complex and then assigned to its corresponding local structure cluster (a class of Betti numbers). This information is fed into the THMM framework for the purpose of classification.

Figure 10. An image pixel is mapped to a cloud of 12 3-D points of the transformed frequency band vector space to build a pixel simplicial complex.

THMM in Change Detection

Change detection is a classification problem with only two classes—Changed and Unchanged—that applies to a data sequence extracted from the pixel of the remote satellite image at two different points in time, T1 and T2. The THMM formalism fits seamlessly with the change detection application at hand. Because a pixel is an input, a VR simplicial complex is extracted from this pixel.

Figure 11. THMM associated with the same pixel of an image at two different points in time, T1 and T2.

he entire sequence is O=O1O2, whereby each Oi=oi1oi2oi12 designates the 12 3-D points associated with a pixel of the image at time Ti(i=1,2). A local latent structure Si mapped to each Oi(i=1,2) corresponds to the Betti number cluster computed from the simplicial complex generated by the VR constructor. In conclusion, the following instantiations are made:

The classes ω1 and ω2 correspond to Changed and Unchanged categories, respectively.

The sequence Oi corresponds to the 12 3-D points generated for a pixel using DWT.

The symbol Si represents the Betti number class (a cluster, aka local structure) assigned to the sequence Oi.A symbol oij of a sequence Oi represents a 3D point among these 12 points.

A hidden state qi is a frequency band Bi selected during the DWT procedure. A set of frequency bands has been selected from which the 12 points have been generated.

Finally, maximum likelihood is applied in the computation of the best ω* among ω1 and ω2 according to the following criterion:


It is worth noting that maximum likelihood-based technique is the one most commonly used for land cover classification especially for medium resolution imagery. However, classification based on maximum likelihood alone has many limitations in separations of remote sensing data such as salt-and-pepper effects and challenges caused by mixed pixels, especially for high-resolution images. To overcome this drawback, maximum likelihood is used only in the computation of the best class ω* among ω1 and ω2 that optimizes the criterion defined in Equation 10. It is worth noting that during the supervised training, the sequences O=O1O2 are tagged changed or unchanged, depending on whether a pixel has changed or not in the image viewed at time T2. The ratio of the number of changed pixels to the total number of pixels in the image represents the magnitude of change. Figure 11 illustrates the instantiation process of the THMM for change detection.


Data Collection

To analyze the efficiency of the proposed model, two change detection data sets were selected, namely the Change Detection dataset and Sztaki AirChange Benchmark. The Change Detection data set comprises 1,000 pairs of images with (800×600)=480,000 pixels each. An image pair corresponds to one reference image and one test (or moving) image. Each pair is associated with a binary ground truth image with the same resolution. The images were constructed by the serious game Virtual Battle Station 2 with Bohemia Interactive Simulations. The data set consists of 100 different scenes containing several thematic (vegetation, water, buildings) and moderate ground reliefs. Each scene was viewed from five different angles, where camera positions were varied at steps of 10 degrees according circle of radius 100 meters at approximately 250 meters high, and with a fixed tilt of about -70 degrees. All images were acquired with a resolution of about 50 cm per pixel. To obtain more information about this data set, please refer to the “Change Detection dataset”. Figure 12 depicts a sample of image pairs at two points in time T1 and T2 and their ground binary truth images, highlighting the changed area with white dots.

Figure 12. Samples of images at time T1 (targets: column 1), same images at time T2 (moving: column 2), and binary ground truth corresponding images (column 3). White dots correspond to the changed areas.

The Sztaki AirChange Benchmark is the second database used. It comprises 13 aerial image pairs of size 952×640 and with a resolution of 1.5 meters and their corresponding ground of truth. Images constituting each pair were acquired at different times over several years, allowing the generation of relevant changes, such as new built-up regions, building operations, planting of large groups of trees, and so on. This database also contains challenging problems in which reversible changes are considered as unchanged, where variations of the vegetation greenery or even shadows are labeled as unchanged in the ground truth. Figure 13 illustrates some samples of image pairs and their corresponding masks (ground truth images).

Figure 13. Samples of images at time T1 (first column), same images at time T2 (second column), and their corresponding ground truth images (third column).

Training Phase

To measure the generalization power of the THMM classifier, a three-fold cross-validation estimation technique is performed in the entire set of 1,000 pairs of pixels. This set is divided into three subsets of 333 pairs of images each, and selected two subsets (666 pairs) for training and one subset for testing (333 pairs). All feature vectors are extracted using the MATLAB platform (MATLAB, 2015). The visible observation sequences (O=O1O2) are extracted from these 666 pairs of images. This procedure is performed using pixels from both classes, Changed and Unchanged, separately. Five Gaussian mixtures (M=5) are used for the estimation of the emission probability values Prob(ot|qt). Therefore, two optimal models 1*=[π1*,A1*,B1*,CC*,D1*] and 2*=[π2*,A2*,B2*,C2*,D2*] are obtained during each training phase. Testing is conducted following each training phase. This procedure is repeated three times with three different subsets for testing. The predicted classification for each pixel is compared to the corresponding pixel in the ground binary truth image. The results in each training/testing procedure are finally averaged to obtain an overall classification performance.

Testing: Decision and Classification

It is well known that a classification problem taking into account a reject option (whereby the classifier has the possibility to hold back its affectation) incurs a cost lower than misclassification (Duda, Hart, & Stork, 2001). Indeed, there are a significant number of cases in which radiometric values of a group of pixels have changed slightly because of a modification in the weather conditions (e.g., presence of clouds). These image regions are often classified as Changed, but no change has actually occurred. The geometric characterization of local structures assigned to these regions can only partially mitigate this false positive effect. A reject threshold can be determined via the risk minimization given by loss matrix. However, the loss matrix element λ representing the loss incurred for choosing the reject option is computed using a selected validation set. Furthermore, when the THMM classifier makes a decision for a test observation sequence, it should take into account both the potential gain and potential loss. An accepted low-risk pixel increases profit, while a rejected high-risk pixel decreases loss. The situation is much more critical and far from symmetric in the area of remote sensing change detection. Denote αi(i=1,2) the action of assigning a single pixel (viewed via O)to class ωi (ω1 = Changed and ω2 = Unchanged, respectively) and α3 the action assigned to the “reject” (or doubt) decision. Likewise, designate λik as the loss incurred for taking action αi when the input actually belongs to class ωk. For an optimal classifier performance, statistical decision theory advocates the action αi with a minimum conditional risk:


A correct classification is assigned a zero loss; however, the classification errors are not equally costly in the remote sensing change detection problem. In fact, wrongly classifying the changed pixels is more critical than misclassifying unchanged pixels. To discriminate wrong decision losses, the parameters λ12 and λ21 are set to 1 and 0.5, respectively. The loss matrix is then expressed as


Therefore, the conditional risk of reject is

{Prob(ωC|O)>2×λ andProb(ωU|O)>λ.

To precisely evaluate the system performance, the receiver operating characteristic (ROC) curve is invoked (Li & Guo, 2014). The results from the folds are then averaged for each performance measure using k-fold cross-validation estimation technique to produce a single estimation.

Experimental Results and Interpretation

The THMM approach has been compared with the traditional HMM. In the HMM framework, a hidden state is assigned a discrete value representing one of the thematic classes (such as “water,” “grass,” etc.). The image was first segmented into thematic classes and each pixel belongs to one thematic class. The THMM has been implemented using the VR approach using the cross-validation criterion that computes the optimal R* with the maximum accuracy (Figure 14). It is worth noting that the classification performance depends on the quality of features and the classifier’s parameters. Different connectivity parameters have been tested during the training phase (parameter tuning phase). The R-value (d=2R) that has produced the highest accuracy (95%) has been selected (R*=0.9).

Figure 14. Cross- validation graph showing the optimal R* value equal to 0.9 that produced the highest accuracy equal to 95% as well as the persistent interval [Birth..Death] = [0.6..0.9] in red.

This selection of a unique d value represents one way to proceed because a filtration (for continuous d-values) of a finite simplicial complex in which persistent homology exhibits qualitative features (barcodes) can also be performed. The choice of selecting a unique d value is justified by the risk of extracting data inherent noise or artifacts. Furthermore, an experimental classification reject threshold value is set to λ=0.02. Figure 15 illustrates an example of change detection result obtained on a test image applied to a change detection data set. Figure 16 depicts a sample of qualitative results over three aerial data sets, namely the Szada, Archive, and Tiszadob. The obtained binary images correspond perfectly to the change mask of the ground truth. Furthermore, these results show that the proposed approach produces homogeneous and smoothing effect. One of the main strengths of the proposed method based on VR simplicial complex representation is its ability to distinguish between real changes and false changes caused by shadows, or seasonal weather effects.

Figure 15. Samples of images at time T1 (targets: column 1), same images at time T2 (moving: column 2), and binary ground truth corresponding images (column 3). White dots correspond to the changed areas.

Figure 16. Qualitative results over three data sets. Column (a), image at T1; column (b), image at T2; column (c), ground truth; column (d), the proposed model results.

For comparison purposes, three classifiers have been implemented: (a) an HMM based only on radiometric values (b) and an HMM that takes into account textural information assigned to a pixel. In order to conduct a comparison between these three classifiers, different values of the classification reject threshold λ that allows different loss matrices to be exhibited are considered. The ROC curves of these classifiers are computed and depicted in Figure 17.

As a second experiment task, the THMM has been compared with other competitive techniques in the field of machine learning (see Dai & Khorram, 1999; Volpi et al., 2013). These techniques were brought inside a MATLAB platform. To ensure reproducibility of the results, this comparative work is conducted on the same database, Change Detection dataset. In other words, some results from the state-of-the-art benchmarked machine learning techniques are provided to derive a global assessment and to widely position the contribution of the THMM formalism. To promote reproducibility and perform a fair model comparison between different change detection approaches, the same set of features prior to classification are extracted.

Figure 17. ROC curves for different values of the reject threshold λ of the HMM (AUC1=0.8548), HMM (considering textural features) (AUC2=0.8748), and THMM (AUC3=0.9083).

Tables 1 and 2 indicate that the THMM outperforms some state-of-the-arts machine learning methodologies when applied to LCCD on both databases (Change Detection dataset and Sztaki AirChange Benchmark). For the purpose of comparison, experiments were conducted on both publicly available data sets, allowing for a fair evaluation. In Figure 17, the area under curve assigned to the THMM is larger than the one assigned to HMMs. The results shown in Tables 1 and 2 clearly underscore the advantage of the THMM over traditional change detection classifiers. The THMM outperforms the other methods because the proposed formalism exploits geometric information, which is highly sensitive to real changes (expressed by geometric representation of pixel sequence). However, these geometric clues cannot always distinguish between real and false changes. According to the geometric formalism, many reversible changes are also labeled as changes. The topological aspect was considered by introducing Betti numbers in extracting features to separate false changes from detected pixels.

Table 1. Comparison of the THMM Approach With Current Machine Learning Techniques on Change Detection Datasets


Change Detector

Overall Accuracy

Dai & Khorram, 1999

Pixel radiometric values at different bands and different times of an image pertaining to the change detection data sets are sequentially fed to the change detector, which is a feed-forward neural network.

Feed-Forward Neural Networks with three layers and four classes (16 change classes).


Volpi et al., 2013

Pixel radiometric values at different bands and different times of an image extracted from the change detection data sets are sequentially fed to the change detector, which is a support vector machine.

The change detection system is a support vector machine with a Gaussian radial basis function as a kernel. The one-against-all scheme is adopted with four classes.



Proposed Methodology

Topological descriptors are embedded within a structural HMM type: They are extracted from the observation sequence using the concept of Betti numbers applied to VR simplicial complexes.

Given two models ω1 and ω2, and a sequence of observations with its topological structure series: compute the model among the two that most likely has generated this sequence.


Note. THMM, Topological Hidden Markov Model; VR, Vietoris-Rips.

Table 2. Comparison of the THMM Approach With Current Machine Learning Techniques on Sztaki AirChange Benchmark


Change Detector


Singh, Kato, Zoltan, and Zerubia, 2014

Information from modified histogram of oriented gradients difference and the gray-level difference

A three-layer Markov random field takes into account information from two different sets of features, namely the modified HOG difference and the gray-level difference.


Benedek and Tamás, 2009

Global intensity statistics with local correlation and contrast features, applying the global energy optimization process to obtain a simultaneously optimal local feature selection

Multilayer conditional mixed Markov model.



(proposed methodology)

Topological descriptors are embedded within a structural HMM type: They are extracted from the observation sequence using the concept of Betti numbers applied to VR simplicial complexes.

Given two models ω1 and ω2, and a sequence of observations with its topological structure series: compute the model among the two that most likely has generated this sequence.


Note. HOG, Histogram of Oriented Gradients; THMM, Topological Hidden Markov Model; VR, Vietoris-Rips.


The THMM presented in this article is a probabilistic land cover change detector that brings about novel information related to the ecological system of a landscape and how people are using this land. Clues about climate change are also inferred from the results obtained. This detector accounts for shapes formed by a sequence of pixels. The topological extension of HMM-based sequence classifiers represents an important step in combining statistics with homology theory in the context of change detection in remote sensing. The experimental results obtained are encouraging because the THMM has outperformed some traditional approaches. This success demonstrates the need for deciphering the conformation of an entire observation sequence. Furthermore, according to the same experiments, the performance of the THMM within the machine learning arena remains significant (around a 4% margin). Future work should focus on determining another criterion for an optimal persistent homology that produces a global (and not local) optimum of the classifier error rate.


We would like to show our gratitude to Dr. Nabil Zerrouki for conducting a part of the experiments in this research.

Further Reading

  • Al-Khudhairy, D., Caravaggi, I., & Giada, S. (2005). Structural damage assessments from Ikonos data using change detection, object-oriented segmentation, and classification techniques. Photogrammetric Engineering and Remote Sensing, 71(7), 825–837.
  • Bouchaffra, D., Cheriet, M., Jodoin, P. M., & Beck, D. (Eds.). (2015). Machine learning and pattern recognition models in change detection [Special Issue]. Pattern Recognition, 48(3), 613–615.
  • Huang, B., Li, Z., & Bo, W. (2009). Spatiotemporal analysis of rural–urban land conversion. International Journal of Geographical Information Science, 23(3), 379–398.
  • Hussain, M., Chen, D., Cheng, A., Wei, H., & Stanley, D. (2013). Change detection from remotely sensed images: From pixel-based to object-based approaches. ISPRS Journal of Photogrammetry and Remote Sensing, 80, 91–106.
  • Inglada, J., & Mercier, G. (2007). A new statistical similarity measure for change detection in multitemporal SAR images and its extension to multiscale change analysis. IEEE Transactions on Geoscience and Remote Sensing, 45(5), 1432–1445.
  • Jianya, G., Haigang, S., Guorui, M., & Qiming, Z. (2008). A review of multi-temporal remote sensing data change detection algorithms. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 37(B7), 757–762.
  • Lu, D., Mausel, P., Brondízio E., & Moran, E. (2004). Change detection techniques. International Journal of Remote Sensing, 25(12), 2365–2401.
  • Ludus, U., Bayarsaikhan, U., Bazartseren, B., Kyung-Ryul, K., Park, K., & Lee, D. (2009). Change detection and classification of land cover at Hustai National Park in Mongolia. International Journal of Applied Earth Observation and Geoinformation, 11(4), 273–280.
  • Masek, J. G., Vermote, E. F., Nazmi, E., Forrest, W., Hall, G., Huemmrich, K.F., . . . Lim, T. (2006). A Landsat surface reflectance dataset for North America, 1990–2000. IEEE Geoscience and Remote Sensing Letters, 3(1), 68–72.
  • Yang, C., & Lo, P. (2002). Using a time series of satellite imagery to detect land use and land cover changes in the Atlanta, Georgia, metropolitan area. International Journal of Remote Sensing, 23(9), 1775–1798.
  • Yuan, F., Kali, E., Loeffelholz, B. C., & Bauer, M. E. (2005). Land cover classification and change analysis of the Twin Cities (Minnesota) metropolitan area by multi-temporal Landsat remote sensing. Remote Sensing of Environment, 98(2), 317–328.
  • Zerrouki N., & Bouchaffra, D. (2014). Pixel-based or object-based: Which approach is more appropriate for remote sensing image classification? In IEEE International Conference on Systems, Man and Cybernetics.


  • Benedek, C., & Tamás, S. (2009). Change detection in optical aerial images by a multilayer conditional mixed Markov model. IEEE Transactions on Geoscience and Remote Sensing, 47(10), 3416–3430.
  • Bengio, Y., & Frasconi, P. (1996). Input-output HMMs for sequence processing. IEEE Transactions on Neural Networks, 7(5), 1231–1249.
  • Bouchaffra, D. (2010). Conformation-based hidden Markov models: Application to human face identification. IEEE Transactions on Neural Networks, 21(4), 597–608.
  • Bouchaffra, D. (2012). Mapping dynamic Bayesian networks to σ-shapes: Application to human faces identification across ages. IEEE Transactions on Neural Networks and Learning Systems, 23(8), 1229–1241.
  • Bouchaffra, D., Govindaraju, V., & Srihari, S. (1999). Post-processing of recognized strings using non-stationary Markovian models. IEEE Transactions on Pattern Analysis and Machine Intelligence, 21(10), 990–999.
  • Dai, X. L., & Khorram, S. (1999). Remotely sensed change detection based on artificial neural networks. Photogrammetric Engineering and Remote Sensing, 65, 1187–1194.
  • De-Alarcon, P. E., Pascual-Montana, A., Gupta, A., & Carazo, J. M. (2002). Modeling shape and topology of low-resolution density maps of biological macromolecules. Biophysical Journal, 83, 619–632.
  • Duda, R. O., Hart, P. E., & Stork, D. G. (2001). Pattern classification. New York, NY: Wiley.
  • Du, P., Xia, J., & Li, F. (2015). Monitoring urban impervious surface area change using China-Brazil earth resources satellites and HJ1 remote sensing images. Journal of Applied Remote Sensing, 9(1), 096094.
  • Edelsbrunner, H. (2014). A short course in computational geometry and topology. Springer Briefs in Applied Sciences and Technology. Cham, Switzerland: Springer.
  • Edelsbrunner, H., & Harer, J. L. (2010). Computational topology: An introduction. Providence, Rhode Island: American Mathematical Society.
  • Fine, S., Singer, Y., & Tishby, N. (1998). The hierarchical hidden Markov model: Analysis and applications. Machine Learning, 32, 41–62.
  • Gong, M., Su, L., Jia, M., & Chen, W. (2014). Fuzzy clustering with modified Markov random field energy function for change detection in synthetic aperture radar images. IEEE Transactions on Fuzzy Systems, 22(1), 98–109.
  • Gong, M., Zhao, J., & Liu, J. (2016). Change detection in synthetic aperture radar images based on deep neural networks. IEEE Transactions on Neural Networks and Learning Systems, 27(1), 125–138.
  • Goodfellow, I., Bengio, Y., & Courteville, A. (2016). Deep learning. Cambridge: MIT press.
  • Gorur, D., & Rasmussen, C. E. (2010). Dirichlet process mixture models: Choice of the base distribution. Computer Science and Technology, 25(4), 615–626.
  • Hajij, M., Wang, B., & Scheidegger, C. (2018). Visual detection of structural changes in time-varying graphs using persistent homology. Proceedings of the IEEE Pacific Visualization Symposium, 125–134.
  • Li, W., & Guo, Q. (2014). A new accuracy assessment method for one-class remote sensing classification. IEEE Transactions on Geoscience and Remote Sensing, 52(8), 4621–4632.
  • Lu, D., Mausel, P., Batistella, M., & Moran, E. (2005). Land‐cover binary change detection methods for use in the moist tropical region of the Amazon: A comparative study. International Journal of Remote Sensing, 26(1), 101–114.
  • Lv, P., Zhong, Y., Zhao, J., Jiao, H., & Zhang, L. (2016). Change detection based on a multifeature probabilistic ensemble conditional random field model for high spatial resolution remote sensing imagery. IEEE Geoscience and Remote Sensing Letters, 13(12), 1965–1969
  • MATLAB. (2015). Version 8.6.0 (R2015b). Natick, MA: The MathWorks Inc.
  • Prasolov, V. V. (2007). Elements of homology theory. Graduate Studies in Mathematics 81, Providence, RI: American Mathematical Society.
  • Rabiner, L., & Juang, B. H. (1993). Fundamentals of speech recognition. Upper Saddle River, NJ: Prentice-Hall.
  • Sanches, I. (2000). Noise-compensated hidden Markov models. IEEE Transactions on Speech and Audio Processing, 8(5), 533–540.
  • Singh, P., Kato, Z., Zoltan, C., & Zerubia, J. (2014). A multilayer Markovian model for change detection in aerial image pairs with large time differences. Proceedings of the IEEE 22nd International Conference on Pattern Recognition (ICPR’14), 924–929.
  • Volpi, M., Tuia, D., Bovolo, F., Kanevski, M., & Bruzzone, L. (2013). Supervised change detection in VHR images using contextual information and support vector machines. International Journal of Applied Earth Observation and Geo-information, 20, 77–85.
  • Zhong, Y., Liu, W., Zhao, J., & Zhang, L. (2015). Change detection based on pulse-coupled neural networks and the NMI feature for high spatial resolution remote sensing imagery. IEEE Geosciences and Remote Sensing Letters, 12(3), 537–541.