- Email: [email protected]
Landsat image interpretation
COMPUTATIONAL STATISTICS & DATAANALYSIS ELSEVIER
Computational Statistics & Data Analysis 20 0995) 75 97
Landsat image interpretation M. Zaki a'*, F. El Tohamy b, A. Hassan c Computer Engineering Dep., AI Azhar University, Egypt b Milita O, Technical College, Egypt Research Center, Egypt
Received June 1993; revised February 1994
The interpretation techniques, used in remote sensing applications, have exploited the wealth of information of the remotely sensed data. These techniques include methods of extracting information from spectral, temporal, and spatial domains of this data. In this paper the remote sensing satellites are considered the primary sources of digital images. Two interpretation algorithms (spectral and temporal) have been designed and implemented on different Landsat images. The first algorithm is directed to the statistical classification of both the Landsat low resolution multispectral scanner (MSS) images and principal component (PC) transformed images whilst the second one is an albedo algorithm to monitor land changes. The algorithm is designed and implemented on MSS and thematic mapper (TM) images. It has been used successfully for land management of the eastern area in Egypt. The results are eventually investigated, argued, and evaluated. Abstract
I. Introduction
Interpretation of digital images Fig. 1 takes different forms in a variety of applications and includes mainly image processing and pattern recognition. Image processing deals with problems in which both input data and output data are pictures. While pictorial pattern recognition deals with methods for producing either a description of the output picture or an assignment of picture to a particular class (Rosenfeld, 1981; Bacher, 198l; Mantas, 1987; Estes and Hajic, 1983; Hall, 1979; Ekstrome, 1984) The largest sources of digital images are the land observational satellites (remote sensing satellites). In fact we are living in the era of satellites since the beginning of *Corresponding author. 0167-9473/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0 1 6 7 - 9 4 7 3 ( 9 4 ) 0 0 0 2 2 - B
76
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
,,.
i in~:U.{J image"1
'ng t
i
p(:
statiztical I classified ,
ol~zitier
I ir~age
restoration eestor.edj " albedo radioM.,qeoM. ie~aqe"I al.ac, rithr,~
land
ohan~es
Fig. 1. Block diagram of image interpretation process.
the 1970s. These digital images were used, through different techniques, to monitor crops and natural resources, prospect for minerals, forecast weather, assess pollution, help planning regions and urban areas, route shipping and in general manage the environment (Estes and Hajic, 1983; Hunt, 1977; Hunt, 1984). Digital image processing includes image restoration, image enhancement, and image compression. While image analysis includes a number of different pattern classification algorithms (Mantas, 1987; Ekstrome, 1984; Woods and Gonzalez, 1981; Fu and Rosenfeld, 1984). A very large number of references are intended to cover different areas of digital image processing field. Rosenfeld (1981) surveyed the image processing and pattern recognition methods and defined the basic differences between them. Mantas (1987) presented an overview of the current methodologies on image pattern recognition. He described the methods of classification processes as a technique of pattern recognition. Many pattern recognition techniques are given in the literature. Fu and Rosenfeld (1984), Horst et al. (1992) and Asrar (1989) have emphasized the main methods of pattern recognition (statistical and structural). Also Fu and Haralick (1983) and Swain (1985) described the main steps of statistical pattern recognition for remotely sensed data. Castleman (1979) and Curran (1988) described the Bayes classification rule. The theorems of image analysis are well-established for mono-chromaticimagery (Llesan*d and Kiefer, 1979). Gonzalez and Wintz (1988) mentioned the classification results of aircraft high resolution MSS and PC transformed images acquired at an altitude of 8000 feet. In this work it will be shown how to make use of the Landsat low resolution data acquired at a mean altitude of 900 km to manage our own earth resources. This work aims at developing an integrated package for land management. The land manager, in question, Fig. 1, is an interpretation-based manager, and it provides solutions for two problems. The first problem is the spectral interpretation of Landsat images. A statistical classifier based on the Bayes rule was designed and implemented on Landsat images of different areas. A comparison with other algorithms has been stated. The second problem is the temporal interpretation of
M. Zaki et al./ Computational Statistics & Data Analysis 20 (1995) 75-97
77
Fig. 2. Landsat MSS5 image of Egypt.
Landsat images acquired at different times. An albedo algorithm was designed and exploited for MSS and TM images of an Egyptian area acquired by Landsat at a time difference of six years. The two procedures use the same type of images (MSS and TM), however, they are complement to each other in the sense that the spectral interpretation gives the space features of an image while the albedo algorithm provides the changes of these features with respect to time. Therefore, the proposed approach could be relied in several application areas when both space and time changes of the image features are needed. Such integration of data analysis has been exploited in this work to find out the vegetation decrease with time at A1-Sharkiya district, Egypt. The Landsat input images include: (1) MSS band 5 and TM band 3 images of A1-Sharkiya, Egypt (Figs. 2 and 3). (2) MSS four band images of Rapid City, USA (Fig. 3). (3) MSS four band images of New York City. USA (Fig. 4). The block diagram of the interpretation processes and the flow diagram of the statistical classification and the albedo algorithms are given. The algorithms are implemented on remote image processing system (RIPS).
2. Spectral interpretation of remotely sensed data One of the widely accepted applications of pattern recognition techniques to image spectral interpretation is the multispectral classification. The process employs multiple images of the same scene, obtained from several spectral regions as
78
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75 97
Fig. 3. Landsat TM3 image of Egypt.
input, and produces as output both displays or maps of scene content on a pixel-topixel basis (Landgrebe, 1981). The basic recognition problem is: given a measurement vector (spectral-bands) X, to determine to which of several classes C~, the measurement belongs. In the statistical classification an optimum decision based on Bayes rule may be made if sufficient information about the point probability function is known (Mantas, 1987; Fu and Rosenfeld, 1984). Often only a set of training samples is given. If these samples are given and labeled according to memberships, this classification is called supervised classification (Curran, 1979; Steiner and Salerno, 1975; Llesand and Kiefer, 1979; Horst et al., 1992). Statistical pattern recognition methods are particularly appropriate for remotely sensed data classification. This is due to the randomness of nature. In such situations statistical pattern recognition methods allow for classifications which are most probably correct (Llesand and Kiefer, 1979). Castleman (1979) provided justification for a histogram or a probability function that can be adequately approximated by a normal (Gaussian) probability density function for one band case the normal density function for class i is given by
1
p(x/Ci)
(2g)1/2o_i exp
-
[-(x- m)~] ~//2
where #i = E[x/Ci] a 2 = E[(x
-
is the mean value of the measurements in class i is the variance of measurements in class i
#i)2/(Ci)]
(1)
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
MSS 1
MSS 3
79
MSS 2
MSS 4 Fig. 4. Four MSS bands of Rapid City.
I n p r a c t i c e ]2i a n d a 2 are u n k n o w n a n d m u s t be e s t i m a t e d f r o m t r a i n i n g samples. F r o m statistical t h e o r y , u n b i a s e d e s t i m a t o r for #i a n d ~r~z are given b y /Ji ~ - -
1 2 )~ x j,
(2)
Fli j = 1
o~ -~ L ~ (~- ~,)~, ni j = l
(3)
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75 97
80
where n~ is the number of training pixels in class i and xj is t h e j t h training pixel for class i. Then the estimated probability function for class i is
~
p(x/Ci) -
exp [
(x -- ~i) 2]j-
(4)
Having made the approximation of normal density function, we need only store the mean and variance for each class in the computer rather than the entire histogram. Conceptually, we can generalize this approach to multispectral data by making use of vector/matrix:
X2 X
i2
=
,
Ui =
n
•
0"i21 ,
~i2n
0"i22
(5)
Zi =
L~tin J
tTinl
tTin2
ffinn
where (X) is the data vector, (Ui) is the mean vector of class i, and (Zi) is the covariance matrix for class i. Then the n-dimensional (multivariate) normal density function can be written from theory of statistics as p(X/Ci)
= (2n),/z
ISi11/2 exp
-- ~ (X -
Iti)r X F I(X
-
I~i)
(6)
where, IS~ I is the determinant of the covariance matrix Zi. S - 1 is the inverse of Si, and ( X - U i ) r is the transpose of the vector (X - Ui). The mean vector and covariance matrix for each class are estimated from training data. Let/~i and ~ are given by fii ~- -
x j,
j=
1, 2 , . . . , n i ,
(7)
?li j = l
_1
- A)(xj
-
(8)
ni j = l
where: ni is the number of training pixels in class i. Two points should be observed in using the assumption of normal density function. First, adequate training samples must be available to estimate the mean and covariance. For n bands the minimum numbers of training pixels are n + 1. For fewer the covariance matrices will be singular, making the computation impossible. The number 100 or higher is better to provide accurate estimates of the class parameters (Haralic and Fu, 1983). Second classes of multimodal histogram cannot be approximated by the normal density function which is unimodal. The solution of this problem is to subdivide the multimodal class into a number of subclasses of unimodal histogram.
M. Zaki et al./ Computational Statistics & Data Analysis 20 (1995) 75-97
81
2.1. The statistical classifier algorithm To classify an image into classes we have to know the probability density of the feature values for each class (from the training sample) and the probability with which the classes occur. Then we can derive a classification criteria that minimize the expected classification error. Let the feature value be denoted by x, and let the probability density of the values of x for class (C, be p(x/Ci) so that k
p(C,) = 1,
i = 1, 2, . . . , k.
(9)
i=1
Then the overall probability density of the values of x for the entire picture is k
p(x) = ~ p(Ci)p(x/C,).
(10)
/=1
So, if we know p(x/Ci) and p(C,), then for any given x we can determine whether it is more likely to arise from a pixel in class C~ or not. Then the joint probability p(C,, x) that a pixel belongs to class C~ and has value x satisfies (from Bayes rule) the following:
p(Ci, x) = p(Ci) p(x/Ci) -- p(Ci/x)p(x),
(11)
where p(x) is the a priori probability that a pixel has value x. It follows that
p(CJx) = p(C,) p (x/C,)/p(x).
(12)
Thus if the p(Ci)'s and p(x/C,)'s are known we can determine the C, for which p(C~/x) is the greatest for any given x. In practice, the p(C,)'s and p(x/C~)'s can be estimated from a given training set for which the pixels have been most probably correctly classified. The classification flow diagram is given in Fig. 7. More generally, in the case of pixel classification on the basis of multifeatures data. Let x be the m-dimensional vector whose components are the given feature values (MSS bands). Then we can define p(X), p(Ci/X), and p ( X / C i ) j u s t as above, and can analogously determine the C, for which p ( C J X ) is the greatest for any given X thus
Vi [ p ( C J X ) = p(C,)p(X/C,)/p(X)].
(13)
The approach discussed in this section minimizes the expected classification error since it assigns each pixel to the most probable class. This approach is applied to remotely sensed data to obtain a pixel by pixel classification. The data represents an area consisting of vegetation, soil, and water covers. The relative area of each cover is calculated using density slice and inventory. Thus the prior class probability p(Ci) = p(vegetation), p(C2) = p(soil) and p(C3) = p(water) are known. The probability p(X) is determined. Then, we calculate p(vegetation/X), p(Soil/X), and p(Water/X). Then X is assigned to the greatest of the above three values.
82
M. Zaki et al./ Computational Statistics & Data Analysis 20 (1995) 75 97
MSS 1
MSS 2
MSS 3
MSS 4
Fig. 5. Four MSS bands of New York. 2.2. Implementation o f the statistical classifier algorithm
Landsat MSS data are ideally suited for the supervised classification because it is quantified, standardized in format and acquired by a stable repetitive platform. The four spectral bands are acquired in registration simultaneously. The four bands are presented to the computer to be classified. An algorithm has been developed to perform the supervised classification techniques. This algorithm was implemented on RIPS. The input data were the multispectral scanner digital image of "New York city" of Fig. 5 acquired in 1981.
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75 97
83
a ~
w m m u
u ~
m l
w
~
m
N
•
: i ' : i i i
i i i i i i ' i i
•
i
.:
..~...~.~ .....
~. . . . . . . . .
~...~ .......................................................
Fig. 6. MSS4 Histogram.
At first a histogram of the image was processed to get an idea about the number of classes in the image. The number of valleys in the histogram determine often, how many classes exist in the image. Fig. 6 shows the histogram of"New York" image of band 1. Three categories have been identified with different discrimination gray level. Then a color image was displayed on RIPS to train the classifier by choosing well identified training set in each category. The RIPS display capability enables us to select, through iterations an optimal training set to improve the classifier accuracy. This set of training pixels constitute less than 1% of the full scene, then the computer can classify automatically the remaining 99% of the scene. The statistical parameters are calculated. The previous training stage produces histograms of the spectral signature of all the pixels, which have been designated in the training set for each class. This
84
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
I
I
START
]
INPUT IMAGE NAME AND NO. OF BAND TO CLASSIFY
DEFINE
/
CLASSES
~
CALL
CLASSIFY
t DISPLAY CLASSIFIED IMAGE
SPECIFY FILES FOR TRAINING
1 DEFI NE TEST PIXELS
T I ~ I NI NG PI XELS
I
1 I
PARAMETERS C A L C . [ FOR TRAIN PI XELS
CALL
I
CLASSIFY
1
EVALUATE ACCURACY
I
F
I
i
CALL HIST. FOR ] TRAINING PI XELS
I
STOP
]
Fig. 7. The classificationflow chart. analysis leads to a successful classification. For example, an abnormally high standard deviation signifies that pixels of different classes have been intermixed in a single class. The classified image is given in Fig. 8.
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75 97
85
Fig. 8. Four band classified image of New York.
By the same way the image of Rapid City was classified using the four MSS band images. Large number of spectral channels require a substantial amount of computer time. But selecting the optimal number of possible combination of spectral channels can be obtained by using the principal component transformation (Jong et al., 1992). The first two principal components from the four components obtained by using PC transform of Rapid City image are classified. The classified output images of the four MSS bands and the two principal components for Rapid City are illustrated in Figs. 9 and 10. From these two figures we conclude that the principle component transformation besides its importance in data compression, it also leads to good interpretational results. The accuracy of the classification of the four MSS bands of Rapid City image is only 2% more than those of classification using two P G images. While acquiring a data reduction of 50% form the original MSS four band images. Figs. 11 and 12 show the classified MS and TM images of Egypt using three and five classes, respectively. The classification algorithm is verified using test data of 500 samples per class. The results are given in Table 1. In this table various images are subject to the classification process. For each image the classification accuracy is given as well as the corresponding number of bands and number of image classes.
86
M. Zaki et al./ Computational Statistics & Data Analysis 20 (1995) 75-97
Fig. 9. Four band classified image of Rapid City.
Fig. 10. Two band classified image of Rapid City.
M. Zaki et al. / Computational Stat&tics & Data Analysis 20 (1995) 75 97
87
b
Fig. 11. Three class classified image of Egypt: (a) MSS (b) TM.
2.3. Results and landsat image classification A statistical classifier has been designed on the basis of Bayes rule (Fu and Rosenfeld, 1984). To apply that statistical classifier on Landsat images we might proceed as follows: (1) display the spectral band image, (2) select training samples using reference data (maps, photos) and zoom for checking,
88
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
b Fig. 12. Five class classified image of Egypt: (a) MSS (b) TM.
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
89
Table 1. The results of Landsat image classification. Classified image
No. of bands
No. of classes
Accuracy (%)
N.y. MSS N.y. MSS N.y. MSS N.y. MSS R.C. MSS R.C. MSS Egypt MSS Egypt MSS Egypt TM Egypt TM
1 2 3 4 4 2 1 1 1 1
3 3 3 3 2 2 3 5 3 5
33.6 56.2 74.7 92.2 96.4 94.5 54.2 50.3 64.5 56.8
(3) subdivide the multimodal classes to unimodal ones, (4) use feature selection to determine the appropriate dimension of data for classification, (5) perform the classifier algorithm, (6) evaluate the classifier accuracy on test samples. The Landsat data classification implemented in the previous section depends mainly on the following factors. (1) The number of spectral bands available for classification. (2) The number of classes to discriminate. (3) The resolution of the sensor factors are evident from the results of Table 1. Although the application of the normal density function to the classification algorithm was subjective, the attained accuracy was highly acceptable because of the learning capability of the developed classifier. The classifier implementation has proved the a priori assumptions of the model.
2.4. Comparison Actually, there are several classification algorithms that can be used, two of them are simple and practically popular. These two algorithms (Curran, 1988) are given below. (1) Parallelepiped algorithm: In which the range of the training set is defined by the highest and lowest pixel value in each band and appears as a rectangular area in the two-feature (band) scatter diagram. Thus, the multidimensional analog of these rectangles yields a Parallelepiped shape (Curran, 1988). (2) Minimum distance to mean algorithm. It comprises three tasks. First, the calculation of the mean pixel value for each class in the training data. Second, the pixels to be classified in the whole data set are assigned the class of their nearest mean vector. Third, if the pixel is farther than an analyst defined distance from any category mean, then it would be classified as unknown.
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
90
Accuracy
(%)
i00
1
5C
1
2
3
Bands
4
Fig. 13. Accuracy comparison of different Classifiers.
The two algorithms have the advantages that they are simple and do not require a priori assumptions for choosing the shape of the density function. However, they suffer from insensitivity to the variance of different categories (Curran, 1988). The performance of the Parallelepiped and minimum distance to mean classifiers has been considered by different authors (Estes and Hajic, 1983; Haralic and Fu, 1983; Llesand and Kiefer, 1979). An accuracy comparison between these two classifiers and the above statistical classifier against the number of multispectral bands of Landsat images is given in Fig. 13. This figure indicates the superiority of the Bayessian algorithm as an accurate classifier. Because of its practical significance, the statistical approach is taken in order to minimize the classification error due to class overlap. In Fig. 13, the accuracy is measured as the number of correctly classified pixels relative to the total number of pixels. On the other hand, the classifiers calculation operations are given as follows:
Statistical: Parallelepiped: Min. distance:
Add/Sub
Mul/Div
Comp
Sqrt
Exp
3 6 6
9 6 6
2 4 2
3
1
3
-
M. Zakietal./ComputationalStat~tics& DataAna~s~20(1995) ~ - 9 7 Time
91
(min)
5
1
3
i- S t a t i s t i c a l
2.5
2- P a r a l l e l e p i p e d 3- M i n i m u m
1
2
3
4
distance
Bands
Fig. 14. Calculation time comparison of different classifiers.
This indicates that the classifiers that did not require a priori assumptions are faster than the statistical classifier, Fig. 14. Therefore, for the applications in which the speed is an important characteristic, the statistical classifier is applied on the image principal components only, in order to speed up the data analysis. This solution provides an effective speed accuracy trade-off.
3. Temporal interpretation of Landsat images The digital images, acquired by remote sensing satellites, for the same area at different times can be likened to multidate (temporal) interpretation. Phenomena, which can be monitored by remote sensors, change over long periods of time, ranging from days to weeks to seasonal or annual time scales. The automated interpretation techniques based on the spectral variation of an area with respect to time, are also being applied in several application areas. Within agricultural context, multidate imagery can help to record the developmental stages of crops. In crop surveys distinct changes during the growing season can permit discrimination on multidate imagery that would be impossible when using any single date image. In addition to increasing separability of category types, change detection can be performed with multitemporal data.
92
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
An accurate and efficient mean of monitoring the changes in the land cover and its renewable resources is vital to updating the factual and practical bases of land management. Field investigation of land changes is expensive due to cost of labor, lodging, and transportation. Computer analysis of satellite images has the advantages of a low unit cost for the data, a reasonable frequency of observation and the ability to cover large areas in sufficient detail for many land management purposes. With the aid of satellite images, some methods may be used to monitor changes in terrain (Llesand and Kiefer, 1979). The simplest method is the visual comparison but it has a high probability of missing subtle and gradual changes. Automated change detection of digital imagery on the other hand may use image differencing, image rationing or albedo approach. The albedo method is applied since it deals with the reflected energy independent on the characteristics of image sensors, the sun elevation angle variation, and the atmospheric haze effect (Landgrebe, 1981; Janza, 1975). 3.1. The albedo concept The albedo method is based on the hypothesis that the percentage of light reflected from the ground (and its cover) will increase if the vegetation cover decreases (in the visible band) and vice versa. (i.e. an increase in density of vegetation will darken the image). The albedo is the percentage of light reflected from the land surface (reflected energy over incident energy) (Janza, 1975). Using sun elevation angle correction and haze removal correction, the value of albedo is a measure of the type of land surface regardless of the time at which the reflectance is measured. Thus albedo changes can be used in temporal interpretation of land changes in a chosen area including increase of vegetation growth, urban areas expansion and desertification problems. Vegetation growth results in a decrease of reflectance, in the visible band, while the arid and urban lands give higher reflectance in the same band (Suits, 1983). 3.2. The albedo algorithm The albedo method, Fig. 15 is applied on two Landsat images of six years difference. Multispectral scanner (MSS) image of May 1978, and thematic mapper (TM) image of July 1984. Both images represent a common area of about 185 km 2 in the eastern part of Egypt. A registration method (Estes and Hajic, 1983; Llesand and Kiefer, 1979) is applied to register the MSS image to the TM image. Fig. 16 shows the MSS, TM, and registered MSS to which the albedo algorithm is applied. In this work a computer program is designed to calculate the albedo change using digital Landsat computer compatible tapes (CCTs) of the previous MSS and TM images of the above described area. Four equations are used to describe the measurement
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
93
TM5 [ IMAGE 1984
IMAGE 1978
l
l
I RADIOMETRIC
I~.DIONETRIC I
CORRECTION
CORRECTI ON
IMAGE I REGISTRATION
ALBEDO
ALBEDO I CALCULATION
CALCULATION
I
]
ALBEDO I CHANGES q i
l,[ NO CHANGEJ ASSIGN LAND ] CHANGE Fig. 15. Flow diagram of albedo algorithm. of reflectivity in a single Landsat band: Bi(av) -- Radi(av). G + O F F ,
(14)
Bi(min) = Radi(min). G + O F F ,
(15)
Radi(av) = RisinO + Hi,
(16)
Radi(min) = Hi,
(17)
where i = 1. . . . . 4; Bi(av) is the average brightness value in the ith band; H~(min) is the minimum brightness value in the ith band; Radi(av) is the average radiance (mw/cm 2 sr) in the ith band as received by the Landsat scanner; Rad~(min) is the minimum radiance in the ith band; G is the gain per unit radiance;
94
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75 97
b
c
Fig. 16. Egypt images (a) MSS (b) TM (c) registered MSS.
O F F is the offset per unit radiance; Ri is the reflectivity in band i (mw/cm 2 sr); 0 is a solar elevation angle as noted in Landsat image annotator; Hi is the atmospheric haze effect in the ith band (mw/cm 2 sr). F r o m the above equations we can derive that Ri =
Bi(av) - Bi(min) G sin 0
(18)
By definition the albedo is the ratio of reflected radiation to incident one. Then Albedo = --Ri= Bi(av) -- Bi(min) Si Si G sin 0 '
(19)
where Si is the Solar radiance in the ith band at the top of the atmosphere. The selected band is the visible one. It corresponds to bands 4 and 5 of MSS and bands 2 and 3 of TM. The radiometric calibrated data of MSS, and T M are given in Table 2 and Table 3, respectively. Both tables report the values of the undertaken
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
95
Table 2. MSS image radiometric calibrated data. L a n d s a t (3): M a y 1978 Band
Wavelength
Lower Reflectance
Upper Reflectance
Gain (G)
Offset
4 5 6 7
0.5-0.6 0.6-0.7 0.7-0.8 0.8-1.1
0 0 0 0
127 127 127 63
58.80 73.84 89.44 14.38
0.04 0.03 0.03 0.03
Table 3. T M image radiometric calibrated data. L a n d s a t (4): July 1984 Band
Wavelength
Lower Reflectance
Upper Reflectance
Gain (G)
Offset
1 2 3 4 5 6 7
0.45 0.52 0.52-0.63 0.63-0.76 0.76-0.90 1.55-1.75 10.4-12.5 2.08-2.35
70 25 20 10 5 123 2
255 193 238 185 255 182 211
230.82 90.67 150.38 80.54 360.62 1 590.94
-0.066 -0.1.57 -0.112 -0.232 --0.086 0 --0.051
Table 4. The results of albedo algorithm. Window No.
Albedo change
Description
1 2 3
no no no
land quality rest as before ,, ,. . . . . . . . . . . . . . . . .
4
no
,,
no no decrease decrease no no no decrease
. . . . . . . . . .
5 6 7 8 9 10 11 12
,,
,,
,,
~,
increase of vegetation growth increase of vegetation growth land quality rest as before . . . . . . . . . . increase of vegetation growth
96
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
band, wavelength and lower and upper reflectance. Also, the gain and effect per unit radiance are given. Table 2 has been recorded on May 1978 while the values of Table 3 are obtained on July 1984. The underlined area was divided into 12 windows (sub-scenes) as shown in Fig. 16 and the albedo value was calculated for each one in both images. The albedo change was calculated. The results are tabulated out from the albedo program. Table 4 reports the achieved results. For each windows the albedo changes indicate whether vegetation or desertification has occurred. No albedo change means that the land quality is not affected.
4. Conclusion
Various algorithms have been presented to interpret effectively Landsat data. In this paper the following techniques are relied upon. (1) A statistical classifier, which is based on Bayes rule, has been used to enhance the classification precision. The accuracy accomplished using MSS images reaches 96%. This value is higher than the accuracy given in all earlier works using minimum distance to mean and paralMepiped classifiers. This achievement is due to (i) The adequate choice of class training samples using the image display capabilities of the RIPS system. (ii) Using prior information to guide the classification process. The high classification accuracy is actually important for interpretation of low resolution remotely sensed data, acquired from high altitudes by land observation satellites. This aspect is an essential feature of the underlied classifier. Such a feature becomes salient when we compare our results with those reported in (Gonzalez and Wintz (1988)), where the images are taken at an altitude of 3000 ft by an aircraft and the resolution is better than 1 m. In both works the classifier performance has the same order of accuracy although, in this work the Landsat images are obtained at a mean altitude of 900 km and a resolution of 80 m. (2) An albedo algorithm to monitor land changes, has been designed and implemented. This algorithm uses MSS and TM images and exploits their temporal interpretation to manage the land resources. These algorithms have been used to make up a Landsat image interpretation package. Such a package is applied successfully to interpret. Landsat images of different areas. In addition, the proposed albedo algorithm has been used to detect the land changes in A1-Sharkiya state, Egypt. This temporal interpretation of different-generation Landsat images provides a cost-effective solution to the problem of land management especially for the developing countries.
References Arvidson, R.E., The earth observing system of the 1990's, IEEE Proc. (June 1985). Asrar, G., Theory and application of optical remote sensing (Wiley, New York, 1989).
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
97
Bacher, K. E., Design of massively parallel processor, IEEE Trans. comput., 29(9) (September 1981). Bernstein, R., Image geometry and rectification, Manual Remote Sensing, 2 (1983). Billengsley, F., Data processing and reprocessing, Manual remote sensing, 2 (1983). Castleman, K., Digital image processing (Prentice-Hall, Englewood Cliffs, NJ, 1979). Conradson, K. and G. Nilsson, Edge enhancement of Landsat images, Comput. vision Image processing, CVIP--42 (May 1988). Curran, P., Principal of remote sensing (Longman, New York, 1988). Ekstrome, M.P., Digital image processing techniques (Academic press, New York, 1984). El Bata, A., Performance of digital image coding systems, Ph.D. Thesis (MTC, Cairo, 1990). Estes, J. and E. Hajic, Fundamentals of image analysis, Manual remote sensing, 2 (1983). Fu, K. and A. Rosenfeld, Pattern recognition and computer vision, IEEE Trans. comput. (October 1984). Gonzalez, R.C. and P. Wintz, Digital image processing (Addison-Wesley, Reading, MA, 1988). Hall, E.L., Computer image processing and recognition (Academic press, New York, 1979). Haralic, R. and K. Fu, Pattern recognition and classification, Manual remote sensing, 2 (1983). Hay, A., Sampling design to test land use map accuracy, Photo. Eng. Remote Sensing, No. 45 (1979). Hixon, M., Evaluation of different schemes for classification of remotely sensed data, Photo. Eng. remote sensing, no. 46 (1980). Horn, B. and R. Woodman, Destripping satellite images, Proc. Image Understanding Workshop (Cambridge, MA, 1978). Hord, R. M., Digital image processing of remotely sensed data (Academic Press, New York, 1982). Horst, B., W. Schneider and A.J. Inz, Multispectral classification of landsat images using neural networks, IEEE Trans. Geoseience remote sensing (May 1992). Hummel, R. and S. Zucker, Histogram modification techniques, Comput. Graphics Image Processing, 5 (1976). Hunt, B.R., Digital image restoration (Prentice-Hall, Englewood Cliffs, NJ, 1977). Hunt, B., Karhunen-Leove multispectral image restoration, IEEE Trans. ASSP, 32(3) (June 1984). Jain, A.K., Advances in mathematical models for image processing, IEEE Proc., 69(5) (May 1981). Janza, F.J., Interaction mechanisms, Manual Remote Sensing, 1 (1975). Jong, L., Sen and K. Hoppel, Principal components transformation of multifrequency polarimetric SAR imagery, IEEE Trans. Geoscience Remote Sensing (July 1992). Kundo, A. and S. Mitra, Application of mean filtering for removal of noises, 1EEE Trans. ASSP, 32 (June 1984). Landgrebe, D.A., Analysis technology for land remote sensing, IEEE Proc. (May 1981). Llesand, T. and R. Kiefer, Remote sensing and image interpretation (Wiley, New York, 1979). Mantas, J., Methodologies in pattern recognition, Pattern recognition society, 20(1) (1987). Moik, J., Digital processing of remotely sensed images (NASA, Washington, DC, 1980). Norwood, V. and J. Lansing, Electro-optical imaging sensor, Manual Remote Sensing, 2 (1983). Oetz, A., Optical remote sensing of the earth, IEEE Proc., 73(6) (June 1985). Pratt, W.K., Digital image processing (Wiley, New York, 1978). Rosenfeld, A., Pattern recognition, IEEE Proc., 69 (May 1981). Rosenfeld, A., Image analysis, progress, and prospects, Pattern Recognition Society, 17 (1984). Rosenfeld, A. and A. Kak, Digital Picture Processing (Academic press, New York, 1982). Short, N., The Landsat workshop (NASA publication No. 1078, Washington, DC, 1982). Simonett, D., The development and principal of remote sensing, Manual Remote Sensing, 2 (1983). Suits, G., The nature of electromagnetic radiation, Manual Remote Sensing, (1983). Swain, P., Remote sensing the quantative approach (McGraw-Hill, New York, 1978). Swain, P., Advanced interpretation technique for earth data, IEEE Proc. 69(5) (June 1985). Steiner, D. and A. Salerno, Remote sensing data system processing and management, Manual remote sensing, 1 (1975). Woods, R. and R. Gonzalez, Real time digital enhancement, IEEE Proc., 69(5) (May 1981).
Computational Statistics & Data Analysis 20 0995) 75 97
Landsat image interpretation M. Zaki a'*, F. El Tohamy b, A. Hassan c Computer Engineering Dep., AI Azhar University, Egypt b Milita O, Technical College, Egypt Research Center, Egypt
Received June 1993; revised February 1994
The interpretation techniques, used in remote sensing applications, have exploited the wealth of information of the remotely sensed data. These techniques include methods of extracting information from spectral, temporal, and spatial domains of this data. In this paper the remote sensing satellites are considered the primary sources of digital images. Two interpretation algorithms (spectral and temporal) have been designed and implemented on different Landsat images. The first algorithm is directed to the statistical classification of both the Landsat low resolution multispectral scanner (MSS) images and principal component (PC) transformed images whilst the second one is an albedo algorithm to monitor land changes. The algorithm is designed and implemented on MSS and thematic mapper (TM) images. It has been used successfully for land management of the eastern area in Egypt. The results are eventually investigated, argued, and evaluated. Abstract
I. Introduction
Interpretation of digital images Fig. 1 takes different forms in a variety of applications and includes mainly image processing and pattern recognition. Image processing deals with problems in which both input data and output data are pictures. While pictorial pattern recognition deals with methods for producing either a description of the output picture or an assignment of picture to a particular class (Rosenfeld, 1981; Bacher, 198l; Mantas, 1987; Estes and Hajic, 1983; Hall, 1979; Ekstrome, 1984) The largest sources of digital images are the land observational satellites (remote sensing satellites). In fact we are living in the era of satellites since the beginning of *Corresponding author. 0167-9473/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0 1 6 7 - 9 4 7 3 ( 9 4 ) 0 0 0 2 2 - B
76
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
,,.
i in~:U.{J image"1
'ng t
i
p(:
statiztical I classified ,
ol~zitier
I ir~age
restoration eestor.edj " albedo radioM.,qeoM. ie~aqe"I al.ac, rithr,~
land
ohan~es
Fig. 1. Block diagram of image interpretation process.
the 1970s. These digital images were used, through different techniques, to monitor crops and natural resources, prospect for minerals, forecast weather, assess pollution, help planning regions and urban areas, route shipping and in general manage the environment (Estes and Hajic, 1983; Hunt, 1977; Hunt, 1984). Digital image processing includes image restoration, image enhancement, and image compression. While image analysis includes a number of different pattern classification algorithms (Mantas, 1987; Ekstrome, 1984; Woods and Gonzalez, 1981; Fu and Rosenfeld, 1984). A very large number of references are intended to cover different areas of digital image processing field. Rosenfeld (1981) surveyed the image processing and pattern recognition methods and defined the basic differences between them. Mantas (1987) presented an overview of the current methodologies on image pattern recognition. He described the methods of classification processes as a technique of pattern recognition. Many pattern recognition techniques are given in the literature. Fu and Rosenfeld (1984), Horst et al. (1992) and Asrar (1989) have emphasized the main methods of pattern recognition (statistical and structural). Also Fu and Haralick (1983) and Swain (1985) described the main steps of statistical pattern recognition for remotely sensed data. Castleman (1979) and Curran (1988) described the Bayes classification rule. The theorems of image analysis are well-established for mono-chromaticimagery (Llesan*d and Kiefer, 1979). Gonzalez and Wintz (1988) mentioned the classification results of aircraft high resolution MSS and PC transformed images acquired at an altitude of 8000 feet. In this work it will be shown how to make use of the Landsat low resolution data acquired at a mean altitude of 900 km to manage our own earth resources. This work aims at developing an integrated package for land management. The land manager, in question, Fig. 1, is an interpretation-based manager, and it provides solutions for two problems. The first problem is the spectral interpretation of Landsat images. A statistical classifier based on the Bayes rule was designed and implemented on Landsat images of different areas. A comparison with other algorithms has been stated. The second problem is the temporal interpretation of
M. Zaki et al./ Computational Statistics & Data Analysis 20 (1995) 75-97
77
Fig. 2. Landsat MSS5 image of Egypt.
Landsat images acquired at different times. An albedo algorithm was designed and exploited for MSS and TM images of an Egyptian area acquired by Landsat at a time difference of six years. The two procedures use the same type of images (MSS and TM), however, they are complement to each other in the sense that the spectral interpretation gives the space features of an image while the albedo algorithm provides the changes of these features with respect to time. Therefore, the proposed approach could be relied in several application areas when both space and time changes of the image features are needed. Such integration of data analysis has been exploited in this work to find out the vegetation decrease with time at A1-Sharkiya district, Egypt. The Landsat input images include: (1) MSS band 5 and TM band 3 images of A1-Sharkiya, Egypt (Figs. 2 and 3). (2) MSS four band images of Rapid City, USA (Fig. 3). (3) MSS four band images of New York City. USA (Fig. 4). The block diagram of the interpretation processes and the flow diagram of the statistical classification and the albedo algorithms are given. The algorithms are implemented on remote image processing system (RIPS).
2. Spectral interpretation of remotely sensed data One of the widely accepted applications of pattern recognition techniques to image spectral interpretation is the multispectral classification. The process employs multiple images of the same scene, obtained from several spectral regions as
78
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75 97
Fig. 3. Landsat TM3 image of Egypt.
input, and produces as output both displays or maps of scene content on a pixel-topixel basis (Landgrebe, 1981). The basic recognition problem is: given a measurement vector (spectral-bands) X, to determine to which of several classes C~, the measurement belongs. In the statistical classification an optimum decision based on Bayes rule may be made if sufficient information about the point probability function is known (Mantas, 1987; Fu and Rosenfeld, 1984). Often only a set of training samples is given. If these samples are given and labeled according to memberships, this classification is called supervised classification (Curran, 1979; Steiner and Salerno, 1975; Llesand and Kiefer, 1979; Horst et al., 1992). Statistical pattern recognition methods are particularly appropriate for remotely sensed data classification. This is due to the randomness of nature. In such situations statistical pattern recognition methods allow for classifications which are most probably correct (Llesand and Kiefer, 1979). Castleman (1979) provided justification for a histogram or a probability function that can be adequately approximated by a normal (Gaussian) probability density function for one band case the normal density function for class i is given by
1
p(x/Ci)
(2g)1/2o_i exp
-
[-(x- m)~] ~//2
where #i = E[x/Ci] a 2 = E[(x
-
is the mean value of the measurements in class i is the variance of measurements in class i
#i)2/(Ci)]
(1)
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
MSS 1
MSS 3
79
MSS 2
MSS 4 Fig. 4. Four MSS bands of Rapid City.
I n p r a c t i c e ]2i a n d a 2 are u n k n o w n a n d m u s t be e s t i m a t e d f r o m t r a i n i n g samples. F r o m statistical t h e o r y , u n b i a s e d e s t i m a t o r for #i a n d ~r~z are given b y /Ji ~ - -
1 2 )~ x j,
(2)
Fli j = 1
o~ -~ L ~ (~- ~,)~, ni j = l
(3)
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75 97
80
where n~ is the number of training pixels in class i and xj is t h e j t h training pixel for class i. Then the estimated probability function for class i is
~
p(x/Ci) -
exp [
(x -- ~i) 2]j-
(4)
Having made the approximation of normal density function, we need only store the mean and variance for each class in the computer rather than the entire histogram. Conceptually, we can generalize this approach to multispectral data by making use of vector/matrix:
X2 X
i2
=
,
Ui =
n
•
0"i21 ,
~i2n
0"i22
(5)
Zi =
L~tin J
tTinl
tTin2
ffinn
where (X) is the data vector, (Ui) is the mean vector of class i, and (Zi) is the covariance matrix for class i. Then the n-dimensional (multivariate) normal density function can be written from theory of statistics as p(X/Ci)
= (2n),/z
ISi11/2 exp
-- ~ (X -
Iti)r X F I(X
-
I~i)
(6)
where, IS~ I is the determinant of the covariance matrix Zi. S - 1 is the inverse of Si, and ( X - U i ) r is the transpose of the vector (X - Ui). The mean vector and covariance matrix for each class are estimated from training data. Let/~i and ~ are given by fii ~- -
x j,
j=
1, 2 , . . . , n i ,
(7)
?li j = l
_1
- A)(xj
-
(8)
ni j = l
where: ni is the number of training pixels in class i. Two points should be observed in using the assumption of normal density function. First, adequate training samples must be available to estimate the mean and covariance. For n bands the minimum numbers of training pixels are n + 1. For fewer the covariance matrices will be singular, making the computation impossible. The number 100 or higher is better to provide accurate estimates of the class parameters (Haralic and Fu, 1983). Second classes of multimodal histogram cannot be approximated by the normal density function which is unimodal. The solution of this problem is to subdivide the multimodal class into a number of subclasses of unimodal histogram.
M. Zaki et al./ Computational Statistics & Data Analysis 20 (1995) 75-97
81
2.1. The statistical classifier algorithm To classify an image into classes we have to know the probability density of the feature values for each class (from the training sample) and the probability with which the classes occur. Then we can derive a classification criteria that minimize the expected classification error. Let the feature value be denoted by x, and let the probability density of the values of x for class (C, be p(x/Ci) so that k
p(C,) = 1,
i = 1, 2, . . . , k.
(9)
i=1
Then the overall probability density of the values of x for the entire picture is k
p(x) = ~ p(Ci)p(x/C,).
(10)
/=1
So, if we know p(x/Ci) and p(C,), then for any given x we can determine whether it is more likely to arise from a pixel in class C~ or not. Then the joint probability p(C,, x) that a pixel belongs to class C~ and has value x satisfies (from Bayes rule) the following:
p(Ci, x) = p(Ci) p(x/Ci) -- p(Ci/x)p(x),
(11)
where p(x) is the a priori probability that a pixel has value x. It follows that
p(CJx) = p(C,) p (x/C,)/p(x).
(12)
Thus if the p(Ci)'s and p(x/C,)'s are known we can determine the C, for which p(C~/x) is the greatest for any given x. In practice, the p(C,)'s and p(x/C~)'s can be estimated from a given training set for which the pixels have been most probably correctly classified. The classification flow diagram is given in Fig. 7. More generally, in the case of pixel classification on the basis of multifeatures data. Let x be the m-dimensional vector whose components are the given feature values (MSS bands). Then we can define p(X), p(Ci/X), and p ( X / C i ) j u s t as above, and can analogously determine the C, for which p ( C J X ) is the greatest for any given X thus
Vi [ p ( C J X ) = p(C,)p(X/C,)/p(X)].
(13)
The approach discussed in this section minimizes the expected classification error since it assigns each pixel to the most probable class. This approach is applied to remotely sensed data to obtain a pixel by pixel classification. The data represents an area consisting of vegetation, soil, and water covers. The relative area of each cover is calculated using density slice and inventory. Thus the prior class probability p(Ci) = p(vegetation), p(C2) = p(soil) and p(C3) = p(water) are known. The probability p(X) is determined. Then, we calculate p(vegetation/X), p(Soil/X), and p(Water/X). Then X is assigned to the greatest of the above three values.
82
M. Zaki et al./ Computational Statistics & Data Analysis 20 (1995) 75 97
MSS 1
MSS 2
MSS 3
MSS 4
Fig. 5. Four MSS bands of New York. 2.2. Implementation o f the statistical classifier algorithm
Landsat MSS data are ideally suited for the supervised classification because it is quantified, standardized in format and acquired by a stable repetitive platform. The four spectral bands are acquired in registration simultaneously. The four bands are presented to the computer to be classified. An algorithm has been developed to perform the supervised classification techniques. This algorithm was implemented on RIPS. The input data were the multispectral scanner digital image of "New York city" of Fig. 5 acquired in 1981.
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75 97
83
a ~
w m m u
u ~
m l
w
~
m
N
•
: i ' : i i i
i i i i i i ' i i
•
i
.:
..~...~.~ .....
~. . . . . . . . .
~...~ .......................................................
Fig. 6. MSS4 Histogram.
At first a histogram of the image was processed to get an idea about the number of classes in the image. The number of valleys in the histogram determine often, how many classes exist in the image. Fig. 6 shows the histogram of"New York" image of band 1. Three categories have been identified with different discrimination gray level. Then a color image was displayed on RIPS to train the classifier by choosing well identified training set in each category. The RIPS display capability enables us to select, through iterations an optimal training set to improve the classifier accuracy. This set of training pixels constitute less than 1% of the full scene, then the computer can classify automatically the remaining 99% of the scene. The statistical parameters are calculated. The previous training stage produces histograms of the spectral signature of all the pixels, which have been designated in the training set for each class. This
84
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
I
I
START
]
INPUT IMAGE NAME AND NO. OF BAND TO CLASSIFY
DEFINE
/
CLASSES
~
CALL
CLASSIFY
t DISPLAY CLASSIFIED IMAGE
SPECIFY FILES FOR TRAINING
1 DEFI NE TEST PIXELS
T I ~ I NI NG PI XELS
I
1 I
PARAMETERS C A L C . [ FOR TRAIN PI XELS
CALL
I
CLASSIFY
1
EVALUATE ACCURACY
I
F
I
i
CALL HIST. FOR ] TRAINING PI XELS
I
STOP
]
Fig. 7. The classificationflow chart. analysis leads to a successful classification. For example, an abnormally high standard deviation signifies that pixels of different classes have been intermixed in a single class. The classified image is given in Fig. 8.
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75 97
85
Fig. 8. Four band classified image of New York.
By the same way the image of Rapid City was classified using the four MSS band images. Large number of spectral channels require a substantial amount of computer time. But selecting the optimal number of possible combination of spectral channels can be obtained by using the principal component transformation (Jong et al., 1992). The first two principal components from the four components obtained by using PC transform of Rapid City image are classified. The classified output images of the four MSS bands and the two principal components for Rapid City are illustrated in Figs. 9 and 10. From these two figures we conclude that the principle component transformation besides its importance in data compression, it also leads to good interpretational results. The accuracy of the classification of the four MSS bands of Rapid City image is only 2% more than those of classification using two P G images. While acquiring a data reduction of 50% form the original MSS four band images. Figs. 11 and 12 show the classified MS and TM images of Egypt using three and five classes, respectively. The classification algorithm is verified using test data of 500 samples per class. The results are given in Table 1. In this table various images are subject to the classification process. For each image the classification accuracy is given as well as the corresponding number of bands and number of image classes.
86
M. Zaki et al./ Computational Statistics & Data Analysis 20 (1995) 75-97
Fig. 9. Four band classified image of Rapid City.
Fig. 10. Two band classified image of Rapid City.
M. Zaki et al. / Computational Stat&tics & Data Analysis 20 (1995) 75 97
87
b
Fig. 11. Three class classified image of Egypt: (a) MSS (b) TM.
2.3. Results and landsat image classification A statistical classifier has been designed on the basis of Bayes rule (Fu and Rosenfeld, 1984). To apply that statistical classifier on Landsat images we might proceed as follows: (1) display the spectral band image, (2) select training samples using reference data (maps, photos) and zoom for checking,
88
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
b Fig. 12. Five class classified image of Egypt: (a) MSS (b) TM.
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
89
Table 1. The results of Landsat image classification. Classified image
No. of bands
No. of classes
Accuracy (%)
N.y. MSS N.y. MSS N.y. MSS N.y. MSS R.C. MSS R.C. MSS Egypt MSS Egypt MSS Egypt TM Egypt TM
1 2 3 4 4 2 1 1 1 1
3 3 3 3 2 2 3 5 3 5
33.6 56.2 74.7 92.2 96.4 94.5 54.2 50.3 64.5 56.8
(3) subdivide the multimodal classes to unimodal ones, (4) use feature selection to determine the appropriate dimension of data for classification, (5) perform the classifier algorithm, (6) evaluate the classifier accuracy on test samples. The Landsat data classification implemented in the previous section depends mainly on the following factors. (1) The number of spectral bands available for classification. (2) The number of classes to discriminate. (3) The resolution of the sensor factors are evident from the results of Table 1. Although the application of the normal density function to the classification algorithm was subjective, the attained accuracy was highly acceptable because of the learning capability of the developed classifier. The classifier implementation has proved the a priori assumptions of the model.
2.4. Comparison Actually, there are several classification algorithms that can be used, two of them are simple and practically popular. These two algorithms (Curran, 1988) are given below. (1) Parallelepiped algorithm: In which the range of the training set is defined by the highest and lowest pixel value in each band and appears as a rectangular area in the two-feature (band) scatter diagram. Thus, the multidimensional analog of these rectangles yields a Parallelepiped shape (Curran, 1988). (2) Minimum distance to mean algorithm. It comprises three tasks. First, the calculation of the mean pixel value for each class in the training data. Second, the pixels to be classified in the whole data set are assigned the class of their nearest mean vector. Third, if the pixel is farther than an analyst defined distance from any category mean, then it would be classified as unknown.
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
90
Accuracy
(%)
i00
1
5C
1
2
3
Bands
4
Fig. 13. Accuracy comparison of different Classifiers.
The two algorithms have the advantages that they are simple and do not require a priori assumptions for choosing the shape of the density function. However, they suffer from insensitivity to the variance of different categories (Curran, 1988). The performance of the Parallelepiped and minimum distance to mean classifiers has been considered by different authors (Estes and Hajic, 1983; Haralic and Fu, 1983; Llesand and Kiefer, 1979). An accuracy comparison between these two classifiers and the above statistical classifier against the number of multispectral bands of Landsat images is given in Fig. 13. This figure indicates the superiority of the Bayessian algorithm as an accurate classifier. Because of its practical significance, the statistical approach is taken in order to minimize the classification error due to class overlap. In Fig. 13, the accuracy is measured as the number of correctly classified pixels relative to the total number of pixels. On the other hand, the classifiers calculation operations are given as follows:
Statistical: Parallelepiped: Min. distance:
Add/Sub
Mul/Div
Comp
Sqrt
Exp
3 6 6
9 6 6
2 4 2
3
1
3
-
M. Zakietal./ComputationalStat~tics& DataAna~s~20(1995) ~ - 9 7 Time
91
(min)
5
1
3
i- S t a t i s t i c a l
2.5
2- P a r a l l e l e p i p e d 3- M i n i m u m
1
2
3
4
distance
Bands
Fig. 14. Calculation time comparison of different classifiers.
This indicates that the classifiers that did not require a priori assumptions are faster than the statistical classifier, Fig. 14. Therefore, for the applications in which the speed is an important characteristic, the statistical classifier is applied on the image principal components only, in order to speed up the data analysis. This solution provides an effective speed accuracy trade-off.
3. Temporal interpretation of Landsat images The digital images, acquired by remote sensing satellites, for the same area at different times can be likened to multidate (temporal) interpretation. Phenomena, which can be monitored by remote sensors, change over long periods of time, ranging from days to weeks to seasonal or annual time scales. The automated interpretation techniques based on the spectral variation of an area with respect to time, are also being applied in several application areas. Within agricultural context, multidate imagery can help to record the developmental stages of crops. In crop surveys distinct changes during the growing season can permit discrimination on multidate imagery that would be impossible when using any single date image. In addition to increasing separability of category types, change detection can be performed with multitemporal data.
92
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
An accurate and efficient mean of monitoring the changes in the land cover and its renewable resources is vital to updating the factual and practical bases of land management. Field investigation of land changes is expensive due to cost of labor, lodging, and transportation. Computer analysis of satellite images has the advantages of a low unit cost for the data, a reasonable frequency of observation and the ability to cover large areas in sufficient detail for many land management purposes. With the aid of satellite images, some methods may be used to monitor changes in terrain (Llesand and Kiefer, 1979). The simplest method is the visual comparison but it has a high probability of missing subtle and gradual changes. Automated change detection of digital imagery on the other hand may use image differencing, image rationing or albedo approach. The albedo method is applied since it deals with the reflected energy independent on the characteristics of image sensors, the sun elevation angle variation, and the atmospheric haze effect (Landgrebe, 1981; Janza, 1975). 3.1. The albedo concept The albedo method is based on the hypothesis that the percentage of light reflected from the ground (and its cover) will increase if the vegetation cover decreases (in the visible band) and vice versa. (i.e. an increase in density of vegetation will darken the image). The albedo is the percentage of light reflected from the land surface (reflected energy over incident energy) (Janza, 1975). Using sun elevation angle correction and haze removal correction, the value of albedo is a measure of the type of land surface regardless of the time at which the reflectance is measured. Thus albedo changes can be used in temporal interpretation of land changes in a chosen area including increase of vegetation growth, urban areas expansion and desertification problems. Vegetation growth results in a decrease of reflectance, in the visible band, while the arid and urban lands give higher reflectance in the same band (Suits, 1983). 3.2. The albedo algorithm The albedo method, Fig. 15 is applied on two Landsat images of six years difference. Multispectral scanner (MSS) image of May 1978, and thematic mapper (TM) image of July 1984. Both images represent a common area of about 185 km 2 in the eastern part of Egypt. A registration method (Estes and Hajic, 1983; Llesand and Kiefer, 1979) is applied to register the MSS image to the TM image. Fig. 16 shows the MSS, TM, and registered MSS to which the albedo algorithm is applied. In this work a computer program is designed to calculate the albedo change using digital Landsat computer compatible tapes (CCTs) of the previous MSS and TM images of the above described area. Four equations are used to describe the measurement
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
93
TM5 [ IMAGE 1984
IMAGE 1978
l
l
I RADIOMETRIC
I~.DIONETRIC I
CORRECTION
CORRECTI ON
IMAGE I REGISTRATION
ALBEDO
ALBEDO I CALCULATION
CALCULATION
I
]
ALBEDO I CHANGES q i
l,[ NO CHANGEJ ASSIGN LAND ] CHANGE Fig. 15. Flow diagram of albedo algorithm. of reflectivity in a single Landsat band: Bi(av) -- Radi(av). G + O F F ,
(14)
Bi(min) = Radi(min). G + O F F ,
(15)
Radi(av) = RisinO + Hi,
(16)
Radi(min) = Hi,
(17)
where i = 1. . . . . 4; Bi(av) is the average brightness value in the ith band; H~(min) is the minimum brightness value in the ith band; Radi(av) is the average radiance (mw/cm 2 sr) in the ith band as received by the Landsat scanner; Rad~(min) is the minimum radiance in the ith band; G is the gain per unit radiance;
94
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75 97
b
c
Fig. 16. Egypt images (a) MSS (b) TM (c) registered MSS.
O F F is the offset per unit radiance; Ri is the reflectivity in band i (mw/cm 2 sr); 0 is a solar elevation angle as noted in Landsat image annotator; Hi is the atmospheric haze effect in the ith band (mw/cm 2 sr). F r o m the above equations we can derive that Ri =
Bi(av) - Bi(min) G sin 0
(18)
By definition the albedo is the ratio of reflected radiation to incident one. Then Albedo = --Ri= Bi(av) -- Bi(min) Si Si G sin 0 '
(19)
where Si is the Solar radiance in the ith band at the top of the atmosphere. The selected band is the visible one. It corresponds to bands 4 and 5 of MSS and bands 2 and 3 of TM. The radiometric calibrated data of MSS, and T M are given in Table 2 and Table 3, respectively. Both tables report the values of the undertaken
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
95
Table 2. MSS image radiometric calibrated data. L a n d s a t (3): M a y 1978 Band
Wavelength
Lower Reflectance
Upper Reflectance
Gain (G)
Offset
4 5 6 7
0.5-0.6 0.6-0.7 0.7-0.8 0.8-1.1
0 0 0 0
127 127 127 63
58.80 73.84 89.44 14.38
0.04 0.03 0.03 0.03
Table 3. T M image radiometric calibrated data. L a n d s a t (4): July 1984 Band
Wavelength
Lower Reflectance
Upper Reflectance
Gain (G)
Offset
1 2 3 4 5 6 7
0.45 0.52 0.52-0.63 0.63-0.76 0.76-0.90 1.55-1.75 10.4-12.5 2.08-2.35
70 25 20 10 5 123 2
255 193 238 185 255 182 211
230.82 90.67 150.38 80.54 360.62 1 590.94
-0.066 -0.1.57 -0.112 -0.232 --0.086 0 --0.051
Table 4. The results of albedo algorithm. Window No.
Albedo change
Description
1 2 3
no no no
land quality rest as before ,, ,. . . . . . . . . . . . . . . . .
4
no
,,
no no decrease decrease no no no decrease
. . . . . . . . . .
5 6 7 8 9 10 11 12
,,
,,
,,
~,
increase of vegetation growth increase of vegetation growth land quality rest as before . . . . . . . . . . increase of vegetation growth
96
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
band, wavelength and lower and upper reflectance. Also, the gain and effect per unit radiance are given. Table 2 has been recorded on May 1978 while the values of Table 3 are obtained on July 1984. The underlined area was divided into 12 windows (sub-scenes) as shown in Fig. 16 and the albedo value was calculated for each one in both images. The albedo change was calculated. The results are tabulated out from the albedo program. Table 4 reports the achieved results. For each windows the albedo changes indicate whether vegetation or desertification has occurred. No albedo change means that the land quality is not affected.
4. Conclusion
Various algorithms have been presented to interpret effectively Landsat data. In this paper the following techniques are relied upon. (1) A statistical classifier, which is based on Bayes rule, has been used to enhance the classification precision. The accuracy accomplished using MSS images reaches 96%. This value is higher than the accuracy given in all earlier works using minimum distance to mean and paralMepiped classifiers. This achievement is due to (i) The adequate choice of class training samples using the image display capabilities of the RIPS system. (ii) Using prior information to guide the classification process. The high classification accuracy is actually important for interpretation of low resolution remotely sensed data, acquired from high altitudes by land observation satellites. This aspect is an essential feature of the underlied classifier. Such a feature becomes salient when we compare our results with those reported in (Gonzalez and Wintz (1988)), where the images are taken at an altitude of 3000 ft by an aircraft and the resolution is better than 1 m. In both works the classifier performance has the same order of accuracy although, in this work the Landsat images are obtained at a mean altitude of 900 km and a resolution of 80 m. (2) An albedo algorithm to monitor land changes, has been designed and implemented. This algorithm uses MSS and TM images and exploits their temporal interpretation to manage the land resources. These algorithms have been used to make up a Landsat image interpretation package. Such a package is applied successfully to interpret. Landsat images of different areas. In addition, the proposed albedo algorithm has been used to detect the land changes in A1-Sharkiya state, Egypt. This temporal interpretation of different-generation Landsat images provides a cost-effective solution to the problem of land management especially for the developing countries.
References Arvidson, R.E., The earth observing system of the 1990's, IEEE Proc. (June 1985). Asrar, G., Theory and application of optical remote sensing (Wiley, New York, 1989).
M. Zaki et al. / Computational Statistics & Data Analysis 20 (1995) 75-97
97
Bacher, K. E., Design of massively parallel processor, IEEE Trans. comput., 29(9) (September 1981). Bernstein, R., Image geometry and rectification, Manual Remote Sensing, 2 (1983). Billengsley, F., Data processing and reprocessing, Manual remote sensing, 2 (1983). Castleman, K., Digital image processing (Prentice-Hall, Englewood Cliffs, NJ, 1979). Conradson, K. and G. Nilsson, Edge enhancement of Landsat images, Comput. vision Image processing, CVIP--42 (May 1988). Curran, P., Principal of remote sensing (Longman, New York, 1988). Ekstrome, M.P., Digital image processing techniques (Academic press, New York, 1984). El Bata, A., Performance of digital image coding systems, Ph.D. Thesis (MTC, Cairo, 1990). Estes, J. and E. Hajic, Fundamentals of image analysis, Manual remote sensing, 2 (1983). Fu, K. and A. Rosenfeld, Pattern recognition and computer vision, IEEE Trans. comput. (October 1984). Gonzalez, R.C. and P. Wintz, Digital image processing (Addison-Wesley, Reading, MA, 1988). Hall, E.L., Computer image processing and recognition (Academic press, New York, 1979). Haralic, R. and K. Fu, Pattern recognition and classification, Manual remote sensing, 2 (1983). Hay, A., Sampling design to test land use map accuracy, Photo. Eng. Remote Sensing, No. 45 (1979). Hixon, M., Evaluation of different schemes for classification of remotely sensed data, Photo. Eng. remote sensing, no. 46 (1980). Horn, B. and R. Woodman, Destripping satellite images, Proc. Image Understanding Workshop (Cambridge, MA, 1978). Hord, R. M., Digital image processing of remotely sensed data (Academic Press, New York, 1982). Horst, B., W. Schneider and A.J. Inz, Multispectral classification of landsat images using neural networks, IEEE Trans. Geoseience remote sensing (May 1992). Hummel, R. and S. Zucker, Histogram modification techniques, Comput. Graphics Image Processing, 5 (1976). Hunt, B.R., Digital image restoration (Prentice-Hall, Englewood Cliffs, NJ, 1977). Hunt, B., Karhunen-Leove multispectral image restoration, IEEE Trans. ASSP, 32(3) (June 1984). Jain, A.K., Advances in mathematical models for image processing, IEEE Proc., 69(5) (May 1981). Janza, F.J., Interaction mechanisms, Manual Remote Sensing, 1 (1975). Jong, L., Sen and K. Hoppel, Principal components transformation of multifrequency polarimetric SAR imagery, IEEE Trans. Geoscience Remote Sensing (July 1992). Kundo, A. and S. Mitra, Application of mean filtering for removal of noises, 1EEE Trans. ASSP, 32 (June 1984). Landgrebe, D.A., Analysis technology for land remote sensing, IEEE Proc. (May 1981). Llesand, T. and R. Kiefer, Remote sensing and image interpretation (Wiley, New York, 1979). Mantas, J., Methodologies in pattern recognition, Pattern recognition society, 20(1) (1987). Moik, J., Digital processing of remotely sensed images (NASA, Washington, DC, 1980). Norwood, V. and J. Lansing, Electro-optical imaging sensor, Manual Remote Sensing, 2 (1983). Oetz, A., Optical remote sensing of the earth, IEEE Proc., 73(6) (June 1985). Pratt, W.K., Digital image processing (Wiley, New York, 1978). Rosenfeld, A., Pattern recognition, IEEE Proc., 69 (May 1981). Rosenfeld, A., Image analysis, progress, and prospects, Pattern Recognition Society, 17 (1984). Rosenfeld, A. and A. Kak, Digital Picture Processing (Academic press, New York, 1982). Short, N., The Landsat workshop (NASA publication No. 1078, Washington, DC, 1982). Simonett, D., The development and principal of remote sensing, Manual Remote Sensing, 2 (1983). Suits, G., The nature of electromagnetic radiation, Manual Remote Sensing, (1983). Swain, P., Remote sensing the quantative approach (McGraw-Hill, New York, 1978). Swain, P., Advanced interpretation technique for earth data, IEEE Proc. 69(5) (June 1985). Steiner, D. and A. Salerno, Remote sensing data system processing and management, Manual remote sensing, 1 (1975). Woods, R. and R. Gonzalez, Real time digital enhancement, IEEE Proc., 69(5) (May 1981).