Automatic change detection in multitemporal remote sensing images using mathematical operations


One of the main tasks of Earth observation is the monitoring of objects and phenomena on its surface. The regularity of the satellite remote sensing allows for constant monitoring and quick detection and evaluation of changes in natural and transformed landscapes. This information is very useful and is applied in land use management, forest management, emergency response and other industries. Regular surveying has been carried out for several decades. During this time a huge archive of satellite imagery has been accumulated. This is, in a way, a chronicle of the state of the Earth’s landcover. It allows you to detect changes years after they occurred, that is, to conduct a retrospective analysis. Among the modern technologies for the processing of satellite images, the technologies of automated detection and cartography of changes in temporal images take an important place. Such technologies are the main element of remote monitoring systems. With this post, we plan to open a series of publications on remote sensing technologies using satellite imagery.



A satellite image is not only a graphical representation of the earth’s surface but also its quantitative description. Each pixel has a quantitative meaning – the value of reflected solar radiation from the surface. Therefore, satellite images can be processed using mathematical methods. Mathematical operations are one of the simplest technologies for detecting changes. These operations are a special case of mathematical processing of raster in general. And there are no specific features. Subtraction and division are used to identify changes. The operation of subtracting the raster can be described by the following formula:


Here, I1 and I2 represent rasters of two dates (I for “image”), D – difference image, index x and у – pixel coordinates, index b indicates multiband image layer, C is a constant added to avoid negative values, but as negative values are not critical it is seldom used.

From the above formula, we see that the pixel value for one date is subtracted from the pixel value with the same coordinates for another date. As a result, the difference image has the same characteristics as the main raster. If these raster layers are of the same number of bands, the difference image will also be multilayer and with the same number of bands (Figure 1).



Fig. 1. The concept of raster subtraction


The operation of dividing space images can be described by the following formula:


Unlike the subtraction formula, the division result is denoted by the letter Q (from the word “quotient”). The constant in the denominator is used to avoid division by zero.


Rasters with a different number of layers can also be subtracted and divided. But for satellite images, it does not make sense, because the images made by one sensor have the same number of channels. And comparing images of different sensors is not the best idea.

When we perform subtraction or division we need to ensure the precision of these procedures. This means that the detected changes in the brightness of the Earth’s surface must correspond to real changes in it. To do this, it’s best to use images taken with one sensor. Then the brightness differences in the images will be associated with changes in the surface, and not with the difference in the imaging technology. If you need to compare the state of the terrain in different years, it is desirable that the images were taken in the same season. Then we are guaranteed to discover long-term changes, not seasonal dynamics of the area.

After selecting the images for comparison, it is necessary to perform their radiometric calibration and atmospheric correction. This way, the difference in the state of the atmosphere will not interfere with detecting real changes in the area.



We will use ENVI to try locating changes, more precisely deforestation, in images. Figure 2 shows two snapshots of the Landsat 5TM images (track/row number 177/025). They cover the same area of deciduous forest, which is located in the regions of the Volchansk (northern part) and Pechenizky (southern part) districts of the Kharkiv region. This forest stretched along the right bank of the Pechenizha reservoir, from Stary Saltov to Pechenigy.

img2left img2right

Fig. 2. Fragments of the Landsat 5 TM image from 1996 (top) and 2007 (bottom)


The first image was acquired on May 7th, 1996. The second – eleven years later – on May 6th, 2007. That is, both images were obtained during the same season, and the pheno-phases of vegetation should be almost identical. Therefore, where there were no changes in the forest cover, the brightness values of the pixels will differ slightly. Of course, the age of the forest has changed. But the corresponding change in brightness will be relatively small and almost identical for all parts of the forest. And where there have been significant changes in the area, changes in brightness will be great. Such changes in the area may include forest fires, deforestation, drying of forests as a result of the spread of diseases and forest pests.

In Fig. 2 we see that forests mostly look the same in 1996 and 2007 images. But for 2007, there are numerous forest cuttings which are not seen in 1996. Our goal is to find such areas in an automated mode.

We will detect changes in the forest area using the image subtraction operation. The division is performed in a similar fashion in ENVI, so we will not look at it separately. The issue of radiometric calibration and atmospheric correction will also not be considered within the framework of this post.  These are separate and quite large themes.


The algorithm for finding deforested areas is shown in Figure 3. It consists of four steps:

1) Create an averaged (pseudo-panchromatic) image from all bands

2) Subtract multitemporal images

3) Create a classification map by slicing

4) Convert the classification map into vector format


The first step requires an explanation. It is not mandatory. When we subtract from one image from another, we get a different image for each channel. Next, you can separately process each difference image. Or you can define the channel for which the difference is represented by the changes in the area best. But in order to simplify the task, the first step of our algorithm is added – averaged images are created from the bands of the multispectral images. Such images are also called pseudo-panchromatic. As with true panchromatic images, they consist of only one layer, are visualized in a black-and-white range and contain information from a wide range of the spectrum. But panchromatic images really acquired by a sensor sensitive to light in a wide range of the spectrum. And pseudo-panchromatic images are made through the averaging of several images in several narrow spectral bands.



Fig. 3. Data processing scheme


1. Creating pseudo-panchromatic images

1) To create an averaged image, select Statistic → Sum Data Bands in the Toolbox. The Sum Data Parameters window pops out (Figure 4). This is a tool that calculates the statistics of a set of images for each pixel and creates corresponding images – sum, mean, standard deviation, and more.

2) In the upper half of the window there is a list of these parameters (Select Outputs Bands:). In our case, only select Mean.

3) The last step is to choose a way to save the result and click OK.
For all the intermediate data processing steps, it’s best to save the results to the temporary memory (in Output Result to, select Memory) and only the final result – the classification map – needs to be saved to file. Then all unnecessary results will be erased from the temporary memory when the program is turned off.



Fig. 4. Setting the calculation of the average image


The results of the calculation of the averaged (pseudo-panchromatic) images are shown in Figure 5. The cuttings that appeared before 2007 are visible as well as in the multi-band image (Fig. 2).


Fig. 5. Images of mean values for 1996 and 2007


2. Creating the difference image

1) To start creating a different image, select Band RatioBand Math in the ENVI Toolbox. After that, a window of mathematical operations for raster will appear (fig. 6, left).

2) The main element of the Band Math window is Enter an expression, which is in the middle of the window. Type here the mathematical expressions for processing rasters. To get a subtraction raster, enter this expression: b2-b1. Letter b denotes variables (images), and the number denotes the channel number. The order of numbering the channels does not matter.

To divide the raster, the expression is b2 / b1. This is the only difference in subtraction and division operations in ENVI.

3) Once the expression is entered, you need to click the Add to List button. If there are no errors in the syntax, the expression will be added to the list (Previous Band Math Expressions) that is located at the top of the Band Math window.

Below the list of expressions, there are four buttons for managing these expressions. The Save button allows you to save the expression in an exp file. You can open an expression from such a file using the Restore button. The Clear button allows you to clear the whole list of expressions together, and the Delete button allows you to remove a particular selected expression from it.

4) If the expression has been added to the Previous Band Math Expressions list, you can click OK. This will bring up a window (Variables to Bands Pairings) in which you need to match the variables from the expression and the raster (Fig. 6, right).

5) The list of variables used in expression is located at the top of the Variables to Bands Pairings window. Below is a list of available ranges (Available Bands list). To match them, you need to left-click to select the first variable in the list and then the raster corresponding to it.



Fig. 6. The task window of the mathematical expression (on the left) and the window for matching the variables with raster (right)


Note that there are no saving options in Figure 6. This is because not all pairs of the variable raster are matched yet. When this is done, the options to save the output (in temporary memory or to a file) will appear in the window.

As a result, we obtain the image of subtraction and division of the images, they are shown in Figure 7. They show clearly visible forest cuts, that look like black areas. Between 1996 and 2007 there was a decrease in the brightness of the surface on their territory.



Fig. 7. Subtraction image (on the left) and division image (to the right) of the investigated area


3. Quantization of the difference image

1) To begin quantization, select the ClassificationRaster Color Slices command in the Toolbox. After that, a window will appear in which we can select the image to be processed. In our case, that is the difference image.

2) After selecting the image for processing, the quantization procedure settings window will appear – Edit Raster Color Slices (Figure 8.9). The default classes will already be set. You need to delete these the крестики button. After that, the window will look like the one in Figure 8. There will be a blank table on the left – a list of classes. The table has three columns: the color of the class, the lower value of the class (Slice Min), the upper value of the class (Slice Max). On the right side of the table is a histogram of the raster values.



Fig. 8. The setup window for the quantization procedure without any classes set


3) You can then set your classes. They are added by clicking the Плюсик button. In our case, we need to add two classes. The first class – deforested areas, the second class – the forested areas that have not undergone major changes.

4) When classes are added, it is necessary to correct the boundaries of the classes, which are indicated in the table. They are set manually, from the keyboard. The optimal values must be determined by the user. After adding classes and configuring their boundaries, the appearance of the window will look like in Figure 9.



Fig. 9. Settings window with the user-specified classes


5) Next, you can adjust the color of the classes. To do this, you must right-click on the appropriate cell in the Color column. A popup menu will appear in which you need to select the Edit Color command … This command opens up the color adjustment window (Figure 10). It allows you to customize the color in different ways:  select from the list, set with the sliders or specify in the color table. In our example, the deforestation was set to red, and for other parts of the forest to blue.



Fig.10. Color setting window


6) After you click the OK button in the quantization window, the classification map will appear as a result. It will be stored in temporary memory. The slicing result for our example is shown in Figure 11. The deforestation with a total area of 216 hectares, representing 2.3% of the forest area, are found.



Fig. 11. Slicing result


7) To save the classification map to a file, you must right-click on the classification chart in the table of contents (Layer Manager). A popup menu will appear, where you can select Export Color SlicesClass Image … or Export Color SlicesShapefile. The first command saves a raster. The second command stores a vector layer in the shapefile format.



Fig.12. Saving slicing results


At this point the process of change detection is over. Saving the output to a shapefile is required for further processing of results in GIS. The use of GIS allows for the creation of quality terrain change maps. You can also conduct in-depth spatial analysis and describe the patterns of location of the detected changes. This is already a step towards prediction of the development of the territory and the provision of measures to prevent unwanted changes.