Partially Supervised Contextual Classification of Multitemporal Remotely Sensed Images Michaela De Martino, Giorgia Macchiavello, Gabriele Moser, Sebastiano B. Serpico Dept. of Biophysical and Electronic Engineering (DIBE) University of Genoa, Via Opera Pia 11a, I-16145 Genoa (ITALY) e:mail:
[email protected]
Abstract — A key-problem in dealing with multitemporal images of a given geographical area is the identification of the changes occurring between distinct acquisition dates. A complete map of the change typologies can be generated when training data are available for all observation dates, but this completely supervised context involves expensive requirements. On the other hand, a completely unsupervised context does not require any prior information but does not allow an analysis of the different typologies of change, since no class information is available at any observation date. In the present paper a contextual multitemporal classification and change detection algorithm is proposed, which deals with remotely sensed image sequences with ground truth information available only at one reference acquisition date. The method integrates clustering information with a two-stage contextual Markov Random Field (MRF) model for the spatio-temporal correlation associated to the sequence. The algorithm is validated on a multitemporal and multispectral real data set, acquired over an agricultural and urban area, and characterized by a large amount of changes between the observation dates. Keywords: Multitemporal classification, change detection, contextual classification, Markov Random Fields.
I.
INTRODUCTION
In the context of environmental monitoring, multitemporal remote sensing images of the geographical area of interest (i.e., sequences of images acquired on the area at different dates) can be useful, for instance, for land cover evolution analysis, for natural resource exploitation or for natural disaster management. In addition, the development of multitemporal analysis techniques plays a strategic role also in view of recent and future missions designed to provide data with short rivisit times (e.g., 12, 24 or 48 hours). A key-problem in dealing with multitemporal imagery is the identification of the changes occurring between consecutive acquisition dates [1, 2]. In particular, an unsupervised change detection methodology allows the identification of changes without prior information for any observation date, which represents an important advantage of the method. However, in a completely unsupervised context, an analysis of the different typologies of change is not feasible, due to the lack of knowledge about the land cover classes at the observation dates. Such an analysis is naturally accomplished in a completely supervised context, i.e., when ground truth data are available for all dates. However, the ground truth availability for all dates involves high costs and is not even realistic in the case of short rivisit time. An interesting
trade-off is thus represented by a partially supervised approach, which exploits ground truth information only for a subset of the acquisition dates. In [3] we developed a partially supervised multitemporal classification and change detection method, dealing with a couple of remote sensing images of the same region. In the present paper we propose a multitemporal classification scheme which extends the aforementioned method and operates with generic image sequences (not only with couples of images) with ground truth information available only for one acquisition date. The method integrates clustering with the use of contextual information, modelling through Markov Random Fields [4, 5] the correlation among neighboring pixels in each image and among distinct images in the sequence. MRFs represent a powerful mathematical tool [4] in image analysis because of their ability to integrate the spatial and temporal contexts in the scene model by the concept of energy functions [4, 6]. Specifically, two MRF-based classification stages are performed here: in the first stage, a spatio-temporal MRF classifier is applied to use the reference image as a temporal context for each of the remaining images of the sequence, in order to fully exploit the available prior information. In the second stage, an MRF classifier is designed, which involves the whole image sequence, by defining a new energy function describing the spatio-temporal correlation of the sequence itself. A comparison between the resulting contextual classification maps allows to generate partially supervised maps of the changes occurring between the acquisition dates. II.
METHODOLOGY
A. Clustering phase Let {I0, I1, …, IN – 1} be a sequence of remotely sensed georeferenced [7] images acquired over the same geographical region at different dates t0, t1, …, tN – 1. We assume a ground truth map G0 to be available only at time t0 (i.e., t0 is assumed as a “reference date”, and I0 as a “reference image”). We propose to identify first the natural classes present in each image Ii (i = 0, 1, …, N – 1) by applying a clustering algorithm. Specifically, the ISODATA method is adopted [8], because it does not require the number K of clusters as an input parameter, but automatically optimizes its choice in a predefined search range. In particular, at date t0, the minimum allowed value for K is identified with the number C of classes in G0 and the maximum allowed value with 5C, in order to take into account the possible presence of multimodal classes: however, this choice is not critical, since it only represents an
This research was funded in part by the Italian Ministry of Education, University and Research. The support is gratefully acknowledged.
0-7803-7930-6/$17.00 (C) 2003 IEEE
upper bound on the number of clusters. For the same reason, the same search range is adopted also at the other dates. Let Mi be the resulting clustering map at date ti, Ki the corresponding number of clusters and Sir the r-th cluster in Mi (r = 1, 2, …, Ki, i = 0, 1, …, N – 1). As in [3], an analysis of the spatial intersection between each cluster in M0 and the training areas in G0 allows to perform a cluster-to-class assignment operation, thus producing a classification map for the reference image. B. First stage MRF model As a first MRF stage, we apply to each couple of images (I0, Ii) (i = 1, 2, …, N – 1) a modified version of the “mutual” MRF classifier proposed in [5] for completely supervised spatio-temporal contextual classification of a couple of remote sensing images. This method cannot be applied directly in the present partially supervised context, since no prior class information is available at date ti for i = 1, 2, …, N. Thus, we adapt the MRF model developed in [5], employing the clusters in place of the ground classes. For each i = 1, 2, …, N, the adopted MRF model defines, for each pixel P ∈ Ii, a 3 × 3 spatial context SCP ⊂ Ii and a 3 × 3 temporal context TCP ⊂ I0 and assigns to the feature vector xP of P, to its cluster label sP ∈ {Sir: r = 1, 2, …, Ki}, and to its context (SCP, TCP) an “energy” value Ui(P). The classification process is then expressed as the minimization of the energy function. The adopted minimization strategy is the Iterated Conditional Mode (ICM) [4], simple and computationally acceptable iterative algorithm, converging to a local minimum of the energy function. This function is assumed to be the sum of three contributions, corresponding to the three considered information sources (spectral information, spatial information and temporal information), i.e.:
U i ( P) = α xiU xi ( xP , sP ) + α siU si ( sP ,SC P ) + α tiU ti ( sP , TC P ) (1) where: Uxi(xP, sP) is a “spectral energy function”, depending on estimates of the posterior probability P(sP| xP) and/or of the cluster-conditional probability density function (pdf) p(xP| sP); Usi(sP, SCP) is a “spatial energy function”, expressing the number of pixels in SCP with label equal to sP; Uti(sP, TCP) is a “temporal energy function” and depends on the Transition Probability Matrix (TPM) [P(Sir| S0k): r = 1, 2, …, Ki; k = 1, 2, …, K0] between the clusters in Mi and the ones in M0: U ti ( sP , TC P ) = −
¦
P( sP | sQ ) ;
(2)
Q∈TC P
αxi, αsi and αti are suitable weight parameters [5]. The cluster-conditional pdf p(xP| sP) is assumed to be Gaussian, and standard sample-mean and sample-covariance are employed as estimates of its parameters. In fact, ISODATA, iteratively minimizing a “sum-of-squared-error” [7, 8] functional, is likely to provide single-mode clusters, for which, especially when dealing with optical data, a parametric Gaussian model is considered to be acceptable [7]. The TPM between the clusters in M0 and the ones in Mi (i = 1, 2, …, N – 1) is estimated through the automatic method developed in [9], which employs the Expectation-Maximization (EM) algorithm, iterative methodology, converging, under mild assumptions, to a maximum likelihood (ML) estimate of the parameters of a probability distribution [10].
C. Second stage MRF model The second stage MRF model aims at systematically exploiting the contextual information present in the whole multitemporal image sequence, by introducing a new energy function, which models the interaction of each image with the previous and the following ones. For each i = 1, 2, …, N – 2 and for each pixel P ∈ Ii, a spatial context SCP ⊂ Ii is defined as in Sec. II.B, and two distinct temporal 3 × 3 contexts TC+P and TC−P are introduced in image Ii + 1 and in image Ii – 1, respectively. Correspondingly, two temporal energy function contributions U ti+ ( sP , TC+P ) and U ti− ( sP , TC−P ) are defined, which describe through the corresponding TPMs, the temporal interaction of Ii with Ii + 1 and Ii – 1, respectively. The resulting energy function is:
U i ( P) = α xiU xi ( xP , sP ) + α siU si ( sP ,SC P ) + +α ti+U ti+ ( sP , TC+P ) + α ti−U ti− ( sP , TC−P )
,
(3)
where all the energy contributions are defined and computed as in Sec. II.B and αxi, αsi, αti+ and αti− are weight coefficients. III. EXPERIMENTAL RESULTS The proposed algorithm is validated on a multitemporal and multispectral data set, acquired over an agricultural and urban area around Pavia, Italy, in April 1994, October 2000 and June 2001 (Fig. 1), respectively. The April image is a a 6-bands Landsat-5 TM image, whereas the October and the June ones are 8-bands Landsat-7 ETM+ images. Due to seasonal land cover modifications and to the long time between April 1994 and the other acquisition dates, a large amount of changes is present in the considered data set. A ground truth map is assumed to be available only at April 1994 and presents five thematic classes, namely “urban”, “wet soil”, “bare soil”, “water” and “wood”. Ground truth maps for the other dates are employed only as test sets: in particular, the October ground map presents the same classes of the April one, whereas in June there is a further “agricultural” class. In order to assess quantitatively the proposed contextual method, the classification accuracy for each thematic class is computed, as well as the overall accuracy (OA) and the average accuracy (AA). The results provided by ISODATA and by the first and the second stage MRF classifiers are listed in Tables 1-3, respectively.
Figure 1. False color composition of optical bands from the data set employed for experiments: April 1994 (left), October 2000 (center), June 2001 (right).
0-7803-7930-6/$17.00 (C) 2003 IEEE
TABLE I. Image
ISODATA: CLASSIFICATION RESULTS Thematic class 0.00%
—
82.07% 74.79%
October 99.62% 97.08% 95.12% 72.11% 22.77%
—
80.75% 77.34%
urban
97.83% 99.65% 93.60% 82.88%
99.66% 91.17% 75.77% 80.91% 59.59% 53.13% 76.54% 76.70%
TABLE II. Image
wood
AA
agr.
June
bare
OA
wet
April
water
FIRST STAGE MRF CLASSIFIER: CLASSIFICATION RESULTS Thematic class
water
bare
wood
urban
wet
October 99.41% 99.11% 89.43% 65.70% 65.19% June
AA
85.13% 83.77%
SECOND STAGE MRF CLASSIFIER: CLASSIFICATION RESULTS Thematic class
water
bare
wood
urban
wet
October 99.89% 95.20% 92.94% 68.51% 74.64% June
—
OA
99.72% 83.14% 98.18% 95.92% 90.17% 29.53% 80.26% 82.78%
TABLE III. Image
agr.
agr. —
OA
AA
86.80% 86.24%
This is further confirmed by a qualitative visual comparison of the corresponding classification maps, which also stresses a better omogeneity of the contextual classification results. The second stage MRF classifier yields a further slight accuracy improvement. The small entity of this improvement can be interpreted as a consequence of the fact that the temporal correlation between couples of images is already exploited in the first stage MRF model and the considered sequence is build up of only three images: in a longer sequence, this contextual classification stage could exploit a larger amount of temporal information, thus more effectively improving the classification accuracy. Finally, the proposed two-stage contextual method sharply improves (with respect to ISODATA) the detection accuracy for the class “wet soil”, thus better separating it from “water”. On the other hand, an accuracy decrease is obtained for “bare soil” and especially for “agricultural” and is interpreted as a consequence of the fact that the values of the parameters of both MRF stages have been manually identified through a “trial and error” procedure. This suggests the importance of the development of automatic parameter setting procedures for such contextual classification methodologies.
99.76% 84.18% 98.79% 97.59% 87.19% 30.67% 80.77% 83.03%
Fig. 2 shows two examples of maps of the typologies of changes between consecutive dates produced by the proposed contextual method. High accuracies are obtained for the classes “water”, “bare soil” and “wood”, whereas worse results are obtained for “wet soil” and “agricultural”, which present a large spectral overlapping with “water” and “wood”, respectively. In particular, it is not possible to distinguish “wet soil” from “water” in April image. An accuracy comparison shows that the proposed first stage MRF classifier yields good accuracies, improving the classification performances with respect to ISODATA. A large increase is obtained especially for the October image (+ 4.38% in OA and + 6.43% in AA). In fact, at this date, ISODATA erroneously labels as “urban” many “wood” and “bare soil” pixels, due to the spectral overlapping between these classes in a purely optical feature space. The first stage MRF model, exploiting the temporal interaction between April and October images, significantly reduces this error.
Figure 2. Examples of change maps generated by the proposed multitemporal second stage MRF classifier: changes from April 1994 to October 2000 (left); changes from October 2000 to June 2001 (right). Color legend: black = no change; blue = “bare soil to wet soil” transition; green = transitions to “agricultural”; white = other changes.
IV. CONCLUSIONS The above experimental results suggest the proposed algorithm to be able to provide quite accurate classification maps for all acquisition dates, despite using training data for only one date , and consequently change maps describing the occurred typologies of changes between consecutive dates. As a future development of this research, we are considering a complete automatization of the parameter setting procedures for the MRF models and an improvement of the clustering results through the adoption of hierarchical clustering procedures. REFERENCES [1]
F. Melgani, G. Moser, S. B. Serpico, "Unsupervised change detection methods for remote sensing images", Optical Engineering, vol. 41, no. 12, pp. 3288-3297, December 2002. [2] A. Singh, "Digital change detection techniques using remotely-sensed data", Int. J. Remote Sensing, vol. 10, n. 6, pp. 989-1003, 1989. [3] G. Moser, F. Melgani, S. B. Serpico, "Partially supervised detection of changes from remote sensing images", Proceedings of IEEE IGARSS2002 Conference, Toronto (Canada), June 2002, vol.I, pp.299-301. [4] R. C. Dubes, A. K. Jain, “Random field models in image analysis”, J. Appl. Statist. (16), 131-163, 1989. [5] F. Melgani, S. B. Serpico, "A mutual approach based on Markov Random Fields for multitemporal contextual classification of remote sensing images", Proceedings of the IEEE IGARSS 2001 Conference, Sydney, Australia, July 2001, vol. VII, pp. 2949-2951. [6] B. Jeon and D. A. Landgrebe, "Classification with spatio-temporal interpixel class dependency contexts", IEEE Trans. Geosci. Remote Sensing, vol. 30, pp. 663-672, Jan. 1992. [7] J. A. Richards, X. Jia, Remote sensing digital image analysis, SpringerVerlag, Berlin, 1999. [8] J. T. Tou, R. C. Gonzalez, Pattern Recognition Principles, AddisonWesley Publ. Comp., Massachusetts, 1974. [9] L. Bruzzone, D. F. Prieto, S. B. Serpico, “A neural-statistical approach to multitemporal and multisource remote-sensing image classification”, IEEE Trans. on Geoscience and Remote Sensing, vol. 37, pp. 13501389, 1999. [10] A. P. Dempster, N. M. Laird, and D. B. Rubin, "Maximum likelihood from incomplete data via the EM algorithm", J. Royal Statist. Soc., vol. 19, n. 1, pp. 1-38, 1977.
0-7803-7930-6/$17.00 (C) 2003 IEEE
