Spatial Change Recognition Model Using Image Processing and Fuzzy Inference System to Remote Sensing

Abstract

After the advent of satellites whose job is to image the surface of the earth, a huge database of imaging data of the surface of the earth was made available to researchers in various sciences to exploit a large data set in their field of work, and the subject of remote sensing gradually came to the attention of researchers in various sciences. For example, geography, environmental science, civil engineering, etc., each analyzed the visual data of the earth’s surface from the perspective of their field. According to this research, the issue of spatial change recognition and their location and calculating the percentage of changes at the ground level has been considered, and the model presented is based on machine vision, image processing, and a fuzzy interface system to reveal features. This research is in the category of applied research and finally, an application will be presented that can lead to the development of software such as Google Earth and can be added to that as an option. Another of the advantages of this model is its easy use compared to specialized software such as Arc GIS, and this is the novelty of this research.

1. Introduction

This paper presents a spatial change recognition model using satellite images, image processing, and a fuzzy interface system as a remote sensing application, and it is in the applied research category. Change recognition in natural phenomena is very important for managing and preserving the environment. On the other hand, video data collection satellites have been able to produce large volumes of images from the surface of the earth and provide them to researchers, such as Google Earth. But the images taken by the satellites show the state of the earth at the time of the image capture, and to get more information about whether or not there has been a change in that area, it is necessary to compare exactly the previous images there. So, change recognition in the study area can give an idea of how to manage and control the environment in that area. Also, change recognition in remotely sensed images is an active research area [1]. Remote Sensing on the Earth’s surface and change recognition mean change detection on the Earth’s surface by processing images of the same geographical area acquired at different times. This field can include forest or vegetation change, forest mortality, defoliation, and damage assessment, wetland change, urban expansion, damage assessment, crop monitoring, changes in glacier mass balance, environmental change and deforestation, regeneration, and selective logging [2]. So, briefly; we refer to some past research in this regard. For example: “Land cover change detection using GIS and remote sensing techniques: A Spatio-temporal study on Tanguar Haor, Sunamganj, Bangladesh” [3], that using a classification area model. Or, “Automated unsupervised change detection technique from RGB color image”, that uses coefficients correlation calculation between color signatures of each two associated pixels from two satellite images for the same area [4], The change detection model has to be insensitive to illumination brightness changes [5], and for this reason, in the proposed model, we use high-resolution grayscale images to prevent changes in light and brightness in satellite images in calculating spatial changes. Or, “Change detection in the city of Hilla between 2007 and 2015 using Remote Sensing Techniques”, that using ArcGIS 10.4 software and those processes include, geometric correction, spectral enhancement, image classification, and cartographic output [6], as mentioned, one of the advantages of this model, that is its easy use compared to specialized software such as Arc GIS, and this is the novelty of this research. Or understanding patterns of vegetation change at the Angkor world heritage site by combining remote sensing results with local knowledge based on analysis stages used to extract spectral plots of pixel values in the region of interest [7], or automatic change detection of buildings in an urban environment from very high spatial resolution images using existing geo-database and prior knowledge based on the Image segmentation model that result of the ratio of the detected area was 86–90% [8], or, On a survey of change detection methods based on remote sensing images for multi-source and multi-objective scenarios was expressed, at present, spatial change recognition as a remote sensing application manifests great significance in numerous change detection applications that are shown on the chart below (Figure 1) [9].

Figure 1.

Published literature statistics of urban change detection according to the keywords remote sensing and urban change detection in Web of Science (total of 1283 publications) [9].

Another study entitled: change detection of soil formation rate in space and time based on multi-source data and geospatial analysis techniques, estimate the dissolution rate and soil formation rate in karst areas of China and analyzed their spatial diversity has been done [10], or in another study, change detection techniques based on multispectral images for investigating land cover dynamics using image processing and mining have been investigated [11]. In another study entitled: Change detection techniques for remote sensing applications, about the distribution of change detection methods [12].

Another study using an images enhancement method, improve the accuracy of SAR¹ images to change detection [14]. Or, in another study, based on similarity measurement between heterogeneous images, change detection has been done [15]. Or based on the maximum entropy principle to obtain the final change detection map and compare with the wavelet-based textural features, plain texture difference, image difference, and log-ratio methods change detection has been done [16].

2. Materials and methods

2.1 Materials

In this research, we want to use publicly available data without the need for special software or special knowledge to collect data and test the proposed model with it. Google Earth is one of the most popular and widely available applications. And since the model is universal and the quality of the images is enhanced by using image processing techniques, it is enough to introduce two images with exactly the same spatial characteristics to the model in two different periods or to display and calculate their differences. So, for example, we take pictures of some different places in the world using the time change feature in Google Earth software and test them with the proposed model. A very important point is that the location specifications are exactly the same in capturing images, which Google Earth software has this feature. So that without changing the desired location, it is enough to change the timeline and capture two images of a place at two different times. Here are examples of images Acquisition from Google Earth from a fixed location at two different time intervals (Figures 2–5).

Figure 2.

A: Greenland 1930, B: Greenland 2021, source and specifications of images: NOAA, US Navy, NGA, GEBCO, Landsat, 74,22′50,50″ N 45,02′01,17″ W, elev 2811 m, height: 3601.47 km.

Figure 3.

A: Jumeirah Palm beach 2000, B: Jumeirah Palm beach 2000, source and specifications of images: NOAA, US Navy, NGA, GEBCO, Landsat, 25,06′42,80″ N 55,03′43,93″ E, elev 10 m, height: 38.26 km.

Figure 4.

A: Oroomiye Lake 1984, B: Oroomiye Lake 2017, source and specifications of images: NOAA, US Navy, NGA, GEBCO, Landsat, 37,08′58,85″ N 45,12′31,05″ E, elev 2318 m, height: 157.76 km.

Figure 5.

South Pole 1957, B: South Pole 2022, source and specifications of images: NOAA, US Navy, NGA, GEBCO, Landsat, 83,32′58,95″ S 62,22′15,95″ E, elev 3178 m, height: 7734.93 km.

Similarly, any other location can be imaged and compared at two intervals, and the proposed model has no limitations, it is enough to keep the exact location specifications and the camera does not change in terms of geographical coordinates or height.

2.2 Methods

The basic model we use to change recognition in images has previously been used to detect changes in industrial parts, where local imaging was performed by the camera and was a type of micrography [17]. The input of the model is the images prepared in the previous section, the type of which was described. In Figure 6, the spatial change recognition model for remote sensing is shown.

Figure 6.

Spatial change recognition model to remote sensing.

2.2.1 Description of model components

2.2.1.1 Acquisition of the first and second spatial image

The source for the acquisition of satellite spatial images with the history is the Google Earth application.

2.2.1.2 Prepare spatial images data

We used the MATLAB image processing toolbox to image mining included, convert RGB images to grayscale and converting the intensity image to double (Pre-processing and data preparation stage), and implemented other parts of the model. The reason for converting images from RGB to grayscale is reducing the data from three dimensions to two dimensions while simplifying the problem. So, in MATLAB, there is a function called “rgb2gray” is available to convert RGB images to grayscale images that we used [17].

2.2.1.3 Edge detection from spatial images with different techniques

In this part, three different methods have been used to edge detection. The reason for using these three methods is that each has its strengths and weaknesses, thus presenting different results in edge detection, which in total will lead to the model’s strength in edge recognition.

2.2.1.3.1 Fuzzy interface systems to edge detection

Briefly, the fuzzy conditions help to test the relative values of pixels that can be present in case of presence on an edge. So, the image is said to have an edge if the intensity variation between the adjacent pixels is large. The mask used for scanning the image is shown in Figure 7 [17].

Figure 7.

Define FIS for edge detection (the first and second spatial image). (fuzzy inference system).

$$\mathrm{Gx}=\left[-11\right],\mathrm{Gy}={\mathrm{Gx}}^{\text{'}}$$

The mask is slid over an area of the spatial image and changes that pixel’s value, and then shifts one pixel to the right and continues to the right until it reaches the end of a row. It then starts at the beginning of the next row and the process continues till the whole image is scanned. When this mask is made to slide over the image, the output is generated by the FIS based on the rules and the value of the pixels [17]. In summary, the steps for using a fuzzy inference system are as follows: a) Crisp spatial images for fuzzified into various FS, having conventional crisp membership functions i.e. Black and White. b) Firing strength is calculated using fuzzy t-norms operators. c) Fuzzy rules are fired for each crisp spatial image. d) Aggregate resultant output FS for all fired rules is achieved by using the max operator (s-norm). e) De-fuzzification using the Centroid method. f) The crisp output is the pixel value of the output image i.e. one containing the edges, and black and white regions. g) The first derivative is performed on the image output from FIS after the application of the noise removal algorithm. h) Further refinement is performed by the second derivative and noise removal [17].

2.2.1.3.2 Sobel’s operator to edge detection

The Sobel operator sometimes called the Sobel–Feldman operator or Sobel filter is used in image processing and computer vision, particularly within edge detection algorithms where it creates an image emphasizing edges [18]. Technically, it is a discrete differentiation operator, computing an approximation of the gradient of the image intensity function. At each point in the image, the result of the Sobel-Feldman operator is either the corresponding gradient vector or the norm of this vector. The Sobel-Feldman operator is based on convolving the image with a small, separable, and integer valued filter in the horizontal and vertical directions. The arrangement of pixels is about the pixel [i, j] shown in Table 1. The Sobel’s operator is the magnitude of the gradient computed by:

${\mathrm{S}}_{\mathrm{X}}=\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ 1& 0& 1\end{array}$

${\mathrm{S}}_{\mathrm{y}}=\begin{array}{ccc}1& 2& 1\\ 0& 0& 0\\ -1& -2& -1\end{array}$

Table 1.

Masks used by Sobel’s operator [18].

$$\mathrm{M}\sqrt{{\mathrm{S}}_{\mathrm{X}}^2+{\mathrm{S}}_{\mathrm{Y}}^2}$$

$${\mathrm{S}}_{\mathrm{X}}=\left({\mathrm{a}}_2+{\mathrm{ca}}_3+{\mathrm{a}}_4\right)-\left({\mathrm{a}}_0+{\mathrm{ca}}_1+{\mathrm{a}}_6\right)$$

With the constant $\mathrm{c}=2$.

Like the other gradient operators, Sx and Sy can be implemented using convolution masks:

2.2.1.3.3 Prewitt’s operator to edge detection

Prewitt’s operator uses the same equations as Sobel’s operator, where constant $\mathrm{c}=1$ (Table 2) [19].

${\mathrm{S}}_{\mathrm{X}}=\begin{array}{ccc}-1& 0& 1\\ -1& 0& 1\\ -1& 0& 1\end{array}$

${\mathrm{S}}_{\mathrm{y}}=\begin{array}{ccc}1& 1& 1\\ 0& 0& 0\\ -1& -1& -1\end{array}$

Table 2.

Masks used by Prewitt gradient operator.

2.2.1.4 Images comparison with different techniques and error reduction

2.2.1.4.1 The structural similarity index measure

The SSIM² formula is based on three comparison measurements between the samples, namely the luminance term, the contrast term, and the structural term. The overall index is a multiplicative combination of the three terms [20].

$$\begin{array}{l}\text{SSIM}\left(\mathrm{x},\mathrm{y}\right)={\left[\mathrm{l}\left(\mathrm{x},\mathrm{y}\right)\right]}^{\mathrm{\alpha }}\cdot {\left[c\left(x,y\right)\right]}^{\beta }\cdot {\left[s\left(x,y\right)\right]}^{\gamma },\mathrm{l}\left(\mathrm{x},\mathrm{y}\right)=\frac{2\text{μxμy}+\mathrm{C}1}{\mathrm{\mu }2\mathrm{x}+\mathrm{\mu }2\mathrm{y}+\mathrm{C}1},\\ \mathrm{c}\left(\mathrm{x},\mathrm{y}\right)=\frac{2\sigma \mathrm{x}\sigma \mathrm{y}+\mathrm{C}2}{\sigma 2\mathrm{x}+\sigma 2\mathrm{y}+\mathrm{C}2}\end{array}$$

$$\mathrm{s}\left(\mathrm{x},\mathrm{y}\right)=\frac{\sigma \mathrm{xy}+\mathrm{C}3}{\sigma \mathrm{x}\sigma \mathrm{y}+\mathrm{C}3}$$

Where $\mathrm{\mu x},\mathrm{\mu y},\sigma \mathrm{x},\sigma \mathrm{y},\text{and}\sigma \mathrm{xy}$ are the local means, standard deviations, and cross-covariance for images x, y. If $\mathrm{\alpha }=\mathrm{\beta }=\mathrm{\gamma }=1$ (the default for Exponents), and $\mathrm{C}3=\mathrm{C}2/2$ (default selection of C3) the index simplifies to:

$$\text{SSIM}\left(\mathrm{x},\mathrm{y}\right)=\frac{\left(2\mathit{\text{μxμy}}+C1\right)\left(2\mathit{\sigma xy}+C2\right)}{\left(\mu 2x+\mu 2y+C1\right)\left(\sigma 2x+\sigma 2y+C2\right)}$$

And, dissimilarity Structural is: (1- SSIM(x,y)).

SSIM measures the perceptual difference between two similar images. It cannot judge which of the two is better: that must be inferred from knowing which the “original” and which has been subjected to additional is processing such as data compression.

2.2.1.4.2 Spatial image subtracts

Each grayscale image is a matrix with color code (0–256), so we can subtract the two matrices to compare the difference between the two images, here we also used the MATLAB image expression box command, Z = imsubtract (img1, img2).

2.2.1.4.3 Absolute difference between the two spatial images

Another way to compare the differences between two spatial images is defined as the sum of the absolute difference at each pixel. The difference value is defined as

$$\mathrm{D}\left(\mathrm{t}\right)=\sum _{\mathrm{i}=0}^{\mathrm{M}}\left|{\mathrm{I}}_{\left(\mathrm{t}-\mathrm{T}\right)}\left(\mathrm{i}\right)-{\mathrm{I}}_{\mathrm{t}}\left(\mathrm{i}\right)\right|$$

Where “M” is the resolution or number of pixels in the image. This method for image difference is noisy and extremely sensitive to camera motion and image degradation. When applied to sub-regions of the image, D (t) is less noisy and may be used as a more reliable parameter for image difference.

$$\mathrm{Ds}\left(\mathrm{t}\right)=\sum _{\mathrm{j}=\mathrm{s}}^{\mathrm{H}/\mathrm{n}}\sum _{\mathrm{i}=\mathrm{s}}^{\mathrm{w}/\mathrm{n}}\left|{\mathrm{I}}_{\left(\mathrm{t}-\mathrm{T}\right)}\left(\mathrm{i},\mathrm{j}\right)-{\mathrm{I}}_{\mathrm{t}}(\mathrm{i},\mathrm{j})\right|$$

Ds (t) is the sum of the absolute difference in a sub-region of the image, where S represents the starting position for a particular region, and n represents the number of sub-regions [21]. There is a function to the MATLAB image processing toolbox to compare two images.

$$\mathrm{Z}=\text{imabsdiff}\left(\mathrm{img}\mathrm{a},\mathrm{img}\mathrm{b}\right)$$

2.2.1.4.4 Histogram comparison

In this part of the model, the histogram of two spatial images is drawn and compared in a map that compares the color changes between 0 and 256 from black to white in spatial images.

2.2.1.4.5 Error calculation and reduction

To check the error rate in the model, it is necessary to compare two images taken exactly the same in terms of space and time from Google Earth with the model. Therefore, ideally, the model should not show any difference. In other words, since the two images taken are exactly the same, the difference calculated by the model must be absolute zero. The accuracy of error calculation is up to 4 decimal places. A comparison of two exactly identical images is shown below (Figure 8):

Figure 8.

Calculate the error rate with an accuracy of 4 decimal places, image description according to Table 3.

2.2.1.5 Information integration

This part, same a single system that information must be shared across all of the functional areas, and the information collected are integrated [22].

2.2.1.6 Spatial changes recognition

We use the following method to calculate the difference between the two spatial images.

$$\text{Percentage difference}=\left(1‐\text{SSIM}\right)\ast 100–\text{Error}$$

And also, we will provide a complete map to show the spatial changes along with calculating the percentage of changes.

3. Result

In this section, after introducing the materials and methods, we run the proposed model to obtain the results. Description of the comparative general map in results is shown in Table 3.

Image number	Description
1	spatial image time1, RGB to grayscale
2	Actions: Fuzzy interface system on the spatial image time1
3	Actions: Prewitt’s operator on the spatial image time1
4	Actions: Sobel’s operator on the spatial image time1
5	spatial image time2, RGB to grayscale
6	Actions: Fuzzy Interface system on spatial image time2
7	Actions: Prewitt’s operator on spatial image time2
8	Actions: Sobel’s operator on spatial image time2
9	spatial image time1 (.jpg) format
10	Difference1, between 6&2, show percentage difference: (1-SSIM) *100
11	Defference2, between 7,3 - apply operator imsubtract
12	Defference3, between 8,4- apply operator imabsdiff

Table 3.

Description of the comparative general map in results.

Experiment 1:

See Figures 9–13.

Figure 9.

A: Greenland 1930, B: Greenland2021, input spatial image to model.

Figure 10.

Comparative general map, the description of the components is as shown in Table 3.

Figure 11.

Comparison of histograms of two temporal spatial images.

Figure 12.

a: Bar histogram of spatial image in time 1, b: Bar histogram of spatial image in time 2, c: Histogram comparison.

Figure 13.

Histogram compression, above axis: The spatial image in time 2 and below: The spatial image in time 1.

Experiment 2:

See Figures 14–18.

Figure 14.

a: Jumeirah Palm beach 2000, b: Jumeirah Palm beach 2000.

Figure 15.

Comparative general map, the description of the components is as shown in Table 3.

Figure 16.

Comparison of histograms of two temporal spatial images.

Figure 17.

a: Bar histogram of spatial image in time 1, b: Bar histogram of spatial image in time 2, c: Histogram comparison.

Figure 18.

Histogram compression, above axis: the spatial image in time 2 and below: the spatial image in time 1.

Experiment 3:

See Figures 19–23.

Figure 19.

a: Oroomiye Lake1984, b: Oroomiye Lake 2017.

Figure 20.

Comparative general map, the description of the components is as shown in Table 3.

Figure 21.

Comparison of histograms of two temporal spatial images.

Figure 22.

a: Bar histogram of spatial image in time 1, b: Bar histogram of spatial image in time 2, c: Histogram comparison.

Figure 23.

Histogram compression, above axis: the spatial image in time 2 and below: the spatial image in time 1.

Experiment 4:

See Figures 24–28.

Figure 24.

a: South Pole 1957, b: South Pole 2022.

Figure 25.

Comparative general map, the description of the components is as shown in Table 3.

Figure 26.

Comparison of histograms of two temporal spatial images.

Figure 27.

a: Bar histogram of spatial image in time 1, b: Bar histogram of spatial image in time 2, c: Histogram comparison.

Figure 28.

Histogram compression, above axis: the spatial image in time 2 and below: the spatial image in time 1.

4. Discussion

Spatial analysis in GIS knowledge is done for experts in this field through specialized software such as ArcGIS, etc., and the issue of recognition of environmental changes, etc. is of special importance to them. In the preparation of image data according to the rules of data mining, the only change in the captured images should be time and other variables of the desired location, including geographical characteristics, camera height from the ground, etc. must remain constant, and otherwise, a calculation error will occur. Therefore, by strictly observing this important point, the system error will be zero; in other words, we will have an accuracy of %100. Also, the test results will be the same no matter how many (repetitions) are performed, which explains the validation of the model.

Comparison segmentation method and proposed model, advantages and disadvantages:

Using segmentation in spatial images is important task for detecting changes. Segmentation must not allow regions of the image to overlap. Thresholding is one of the oldest methods used for image segmentation. It is based on the gray level intensity value of pixels. The histogram of an image consists of Conceptually is similar to the classifiers except that they are implemented in the spatial domain of an image rather than in a feature space. It treats the segmentation as a registration process. Some researchers used atlases not only to impose spatial constraints but also to provide probabilistic information about the tissue model. The advantage is that it can segment an image with no well-defined relation between regions and pixels. K-means is a clustering method that partitions the n-points into the k-clusters in which each pixel belongs to one cluster by minimizing an objective function in such a way that within a cluster sum of squares is get minimized. It starts with k-clusters and each pixel is assigned to one cluster. The limitation of the K-means algorithm is computational time increases on implementation in large amounts of data but our proposed model is independent of clustering so will be faster than the K-means algorithm.

5. Implication

In the proposed model, using machine vision, image data processing, fuzzy mathematical techniques, the use of known masks, as well as historical images in Google Earth software, we can measure changes in images and measure them. This software can be added to Google Earth software as a development part and users can easily view spatial changes in terms of time changes.

6. Conclusion

The change recognition model presented by the author of the article earlier [17], was also used in spatial change recognition. The most important difference between the previous models in change recognition in industrial parts with the current model is in the field of capturing images and image preparation. Here the images are taken from Google Earth and no filter is used to prepare them because the images are of sufficient quality and adding any filter will cause a computational error in the model (this was tested many times by the authors) but in the base model, the capture of images was with a local camera that had the error and the macrographic imaging technique was done [17].