An Incongruence-Based Anomaly Detection Strategy for Analyzing Water Pollution in Images from Remote Sensing

Abstract

The potential applications of computational tools, such as anomaly detection and incongruence, for analyzing data attract much attention from the scientific research community. However, there remains a need for more studies to determine how anomaly detection and incongruence applied to analyze data of static images from remote sensing will assist in detecting water pollution. In this study, an incongruence-based anomaly detection strategy for analyzing water pollution in images from remote sensing is presented. Our strategy semi-automatically detects occurrences of one type of anomaly based on the divergence between two image classifications (contextual and non-contextual). The results indicate that our strategy accurately analyzes the majority of images. Incongruence as a strategy for detecting anomalies in real-application (non-synthetic) data found in images from remote sensing is relevant for recognizing crude oil close to open water bodies or water pollution caused by the presence of brown mud in large rivers. It can also assist surveillance systems by detecting environmental disasters or performing mappings.

1. Introduction

The potential applications of computational tools for analyzing remote sensing data taken from environments containing large rivers attract much attention in the scientific literature [1,2,3,4,5,6,7,8,9,10,11,12]. Two of these tools, anomaly detection [13,14] and incongruence [15,16,17,18], are important for studies to find and categorize non-conforming patterns present in images, such as outliers [19,20,21,22,23,24]. For example, Qi Liu et al. in [25] presented a framework to detect outliers and anomalous events [19,20,21,22,23,24]. They highlighted the importance of basing anomaly detection studies on taxonomies of anomalies. Kittler et al. in [13] proposed a taxonomy of anomalies which expanded the concept of anomaly beyond the conventional meaning of outlier. They used sensory data quality assessment [26], contextual [9] and non-contextual [27,28,29,30] classifiers [5,15,16,17,18], and an incongruence indicator [16,17,18] to identify each type of anomaly [13]. According to this taxonomy [13] anomalies can be, for example, of the types unknown object, measurement model drift, unknown structure, unexpected structural component, component model drift, and unexpected structure and structural components. This taxonomy [13] is well-known and widely accepted by the scientific community because it has the potential to be applied for solving problems in many different research areas [13]. Therefore, studies related to the application of the taxonomy [13] on synthetic data are common, such as in [16,17,18]. However, to the knowledge of the authors, studies have not addressed the practical application of the taxonomy [13] to solve real-world problems [2,31,32,33], because it remains a challenge for all research areas.

Anomaly detection [13] and incongruence [15,16,17,18] are two powerful computational tools from pattern recognition (PR) [3,11,12,13,14,15,16,17,18,34,35] and computer vision (CV) [3,36]. PR is a scientific area of study dedicated to analyze patterns and regularities in data [3]. PR provides powerful tools [34] for many different applications and research areas [35] such as scientific research, private and public industries, military activities, etc., [14,21,22,23,24,32,37,38,39,40,41,42,43,44,45,46,47,48]. For example, PR is important for geosciences as its tools are used to analyze geographical features of environments in digital images from remote sensing, i.e., scenes [11,12,49,50,51,52,53,54]. Additionally, PR also provides powerful tools to help machine perception for CV [3,36]. CV is a scientific area of study committed to developing artificial systems responsible for obtaining information from multidimensional data, such as images [36]. A relevant example of machine perception is the ability of computers to recognize behaviors of geographical features of interest, such as rivers, lakes, fields, plantations, forests, etc. Our interest goes further than a solution based on PR and CV for a human user of a computer to apply in order to solve a specific and already-known problem from remote sensing, such as only detecting water pollution [55,56,57,58,59,60] as an outlier in an image. In this case, we could use spectral analysis [61] (spectral index [62] or slope ratio [63]) to solve our problem directly. Our motivation has been to find a real application to the taxonomy [13] inspired by the real-world situations captured in images from remote sensing [2,31,32,33].

When analyzing an image for the first time, a machine could apply each of all known solutions already used in remote sensing (e.g. spectral analysis [61]), one at a time, until detecting a problem. However, this could be considered trial and error learning [35]. The use of the taxonomy [13] can bring to the computational systems a more organized and structured approach of machine learning [28,35,64,65,66,67,68,69,70,71,72,73,74] for detecting problems, i.e., more “intelligence.” In this case, for example, a decision-making system should itself first identify a problem with the minimum interaction with humans before offering an adequate solution that could be based on spectral analysis [61]. A decision-making system [13,16,17,18] is a software installed in a computer which needs to identify different types of problems in order to automatically select appropriate responses to deal with each type of problem [13]. In other words, the machine should choose the appropriate response to the detected problem based on the type of anomaly identified by the system that could be, for example, unexpected structure and structural components [13]. Our research is applied in this context.

The innovation of this work is to introduce the first solution to a practical application of the taxonomy [13] to solve real-world problems [2,31,32,33]. What makes our study unique is that our study allows computers to categorize real-world problems based on the use of a thorough taxonomy [13] which is a “scaffolding” [75] to increase the machine’s abilities [35]. In other words, beyond identifying problems, this study opens precedents for the machines to recognize the contexts in which the problems are inserted. Each type of anomaly represents a different context [13]. Therefore, the most appropriate response to a real-world problem depends on the context in which that problem is inserted [13]. Humans can learn how to solve problems inserted in different contexts much better than machines [75]. However, our study focuses on machines rather than humans. Providing conditions for a machine to learn how to distinguish and detect different real-world problems taking into consideration the different contexts and based on a single concept (the taxonomy [13]) make our study unique. This leads to a higher level of meaningful learning, i.e., an evolution in machine learning [28,35,64,65,66,67,68,69,70,71,72,73,74].

Among all relevant areas of research, we have found the first solution to this problem just in remote sensing. It opens up an entirely new perspective and a set of possibilities for future work in remote sensing and many other areas. In remote sensing, for example, the already known solutions for specific problems, such as those based on spectral analysis [61], can be used to spread the taxonomy [13] to be applied on different problems. In this research, the real-world situation captured in images from remote sensing which inspired us to investigate solutions in PR and CV is water pollution, mainly regarding occurrences of disasters [31,32,33]. If this situation could be revealed by the disagreement between classifiers [13,15,16,17,18], then water pollution could be discovered by anomaly detection [13].

Therefore, this paper presents an incongruence-based anomaly detection strategy for analyzing images from remote sensing, with the aim of describing its practical application in detecting water pollution. Therefore, the practical application of the taxonomy [13] is introduced as a strategy for detecting anomalies in real-world application [31,32,33], i.e., non-synthetic data present in images from remote sensing. We used divergence [16,17] between two image classifications [5,13,15,16,17,18] (contextual [9] and non-contextual [8,27,28,29,30]) in order to detect incongruence [13,15,16,17,18], revealing the presence of anomalies [13]. This research focuses only on the detection of anomalies of the type unexpected structure and structural components, which have been described by Kittler et al. in [13].

This paper is organized as follows: Section 2 presents a background; Section 3 comments on related work; Section 4 summarizes the proposed strategy; Section 5 describes the materials and methods used for performing the experiments; Section 6 reviews the achieved results; Section 7 discusses the results and Section 8 presents the conclusions.

2. Background

2.1. The Differences between Outlier, Anomaly, and Incongruence

In a statistical analysis, an outlier [19,20,21,22,23,24] is an observation which deviates significantly from other observations [19]. The deviation makes the observation appear to be generated by a different mechanism. Therefore, the observation appears to be inconsistent with the remainder of the data series. For example, consider an image of a large river for which all pixels are blue except one red pixel. In this case, the red pixel is an outlier. Outlier, in the conventional meaning, is a term used to define any type of anomaly [19]. The term “detection of outlier” [19,20,21,22,23,24] must be used when it is not necessary to categorize the anomaly. The term “detection of anomaly” [13] must be used only when it is necessary to categorize the anomaly. Therefore, it is common to find inadequate exchange of these terms in the scientific literature [14,19].

Anomaly detection is the activity of discovering and categorizing non-conforming patterns and occurrences when there is a failure to relate observed sensory data to information we already know about the image [13]. For example, the detection of absence of water in the part of an image that is correspondent to a large river is an out-of-pattern occurrence and, therefore, an anomaly. Because there are different types of anomalies [13], their detection must distinguish each type of anomaly.

For example, anomalies of the type unexpected structure and structural components are those ones detected when the analyses are performed taking into consideration the features which allow classifiers to identify components of samples and how they are structured to compose the samples [13]. The same analyses are applied to detect where differences related to components or their structures can be found in the image, i.e., spatial analyses [13]. Nonetheless, it does not consider when they occurred, i.e., temporal analyses [13]. Therefore, there is no time-series [76,77,78], i.e., satellite images before and after the occurrence of an event [19] or disaster [31,32,33], involved in the context. Moreover, detection of anomaly of this type depends on other three conditions to be satisfied: high sensory data quality [26], contextual [9], and non-contextual [8,27,28,29,30] classification [5,13,15,16,17,18] of the same image, and incongruence [13,15,16,17,18]. An example of this type of anomaly is related to detecting brown mud where there are high levels of turbidity in large rivers. In contrast, for example, incongruent results from contextual and non-contextual classification of images qualified with low values indicate the presence of another type of anomaly [13]. In another example, incongruent results from contextual and non-contextual classification of time-series images [76,77,78] qualified with high values indicate the presence of a third type of anomaly [13].

Incongruence [13,15,16,17,18] is the presence of discordant results and occurs whenever there is contradictory evidence provided by the analyses of sensory data performed by different types of classifiers [13]. For example, taking into consideration the two parallel analyses of the same part of an image that is correspondent to a large river, there is incongruence if water is detected by one analysis and not detected by the other one.

2.2. Non-Contextual and Contextual Classifiers

An environment can be modeled by some computational tools, e.g., those based on machine perception, using samples (aggregation of pixels) from the application domain and previous information about the problem [15]. Models trained with different samples are able to be represented by classifiers [16]. In PR, classifiers are experts or learners which perform classifications, i.e., they assign objects to one class or category among predefined others [79]. While predicting classes for a common input, such as a scene, classifiers are expected to present similar probability estimates despite the environment [16]. Consequently, computational tools, which base their decision-making actions on the use of multiple classifiers, usually expect all classifiers to support the same hypothesis [17].

There are two different types of classifiers, non-contextual and contextual [13]. Non-contextual classification is a general-level task that results in a weak learner, but it is less constrained than the contextual one [79]. It is more powerful than a random classifier, despite being much less precise when compared to the contextual one [79]. Contextual classification is a specific-level task which leads to a strong learner, but it is more dependent on specific knowledge such as previous knowledge or training data [9]. It is composed of multiple weaker classifiers that work in synergy to reach a more robust classification [9].

On one hand, when either two non-contextual [27,28,29,30] or two contextual [9] classifiers diverge in their classifications [13,15,16,17,18], this condition only exposes the limited classification of one of the classifiers. It does not present any relation to incongruence detection [13,15,16,17,18] for both classifiers. On the other hand, models of an environment which are represented by non-contextual and contextual classifiers can be used for identifying incongruences if they are used concurrently [15]. After receiving input data, both types of classifiers produce class-posterior probabilities as results [15]. Incongruence happens when the input data presented to both classifiers induces a significantly large discrepancy between posterior probabilities, i.e., incongruence is a conflicting prediction that happens when the probability calculated by some non-contextual classifier is much bigger than the probability calculated by some contextual classifier [15].

Therefore, the concurrent use of non-contextual and contextual classifiers allows systems to recognize the occurrence of conflicting classifications [15]. It enhances the capability of decision-making systems to detect incongruence [15,16,17,18]. Consequently, decision-making systems can have their control actions also conditioned by the monitoring and detection of incongruence [15,16,17,18]. Decision-making systems can guide computers to give adequate and fast responses after detecting problems [13], e.g., to prevent water pollution disasters or to analyze their extents in the case of their occurrences. For this reason, incongruence is of considerable interest for decision-making systems, since the strategy reveals the occurrence of some type of anomaly [13].

3. Related Work

3.1. Research on Outliers Detection

Most current research focuses on outlier detection, such as [1,14,19,20,21,22,23,24,25,44,51,52,53,54,66,76,77,78,80,81,82], however we discuss below only those research which are more relevant to this article. An outlier-based change detection method to detect abnormal points in multi-scale time-series images from remote sensing was proposed by Yin Shoujing et al. in [76]. The method analyzes time-series images from remote sensing and detects temporal and spatial changes. For that research, changes caused by weather condition variation, climate changes, sensor aging, vegetation phenology, human activities, and emergencies (such as drought, fire, pest, insect, etc.,) were considered abnormal. The research focused on land cover changes. These changes are relatively rare over large areas, even for a long period. The research concluded that the land cover changes registered in the images were outliers.

An algorithm capable of identifying outliers in the global data (distributed earth science databases) without moving the whole data to a single database was presented by Bhaduri et al. in [80]. The algorithm analyzes the data to search for extremes or outliers. The algorithm detects outliers which are missed when data is available only at a single database. Their achieved results were compared to the achievements of other algorithms which also detect outliers. An unsupervised outlier detection framework was proposed by Qi Liu et al. in [25]. The framework requires no prior knowledge to detect outliers and anomalous events. The research defined outliers as objects which have either low spatial or temporal coherence with their neighbors, and it defined an anomalous event as a group of outliers that share similar spatial and temporal anomalous behaviors. The research has defined four different categories of outliers or anomalous events and has reinforced the importance of basing anomaly detection studies on taxonomies of anomalies. Nonetheless, outlier detection is not enough for systems to categorize many types of problems. Anomaly detection can overcome this drawback.

A new daily snow cover dataset was developed by Bormann et al. in [81]. The dataset is a satellite-based observational record useful to characterize snow duration and cover extent. These snow cover observations are important to detect outliers such as anticipated melt of snow, declines in snow cover extent, and short season duration. The dataset was assessed using snow detection estimates derived from Landsat Thematic Mapper data and was compared to another snow cover dataset. Although the study focused on satellite data for the alpine region in Australia, the approach can be applied to other regions and other sensors to help assess snow monitoring. Consequently, the study contributes to research on water resource management and snow hydrology.

A method for detecting anomaly regions in each image of satellite image time-series was proposed by Zhou et al. in [77]. The method identifies spatial-temporal dynamic processes of unexpected changes of land cover based on seasonal autocorrelation analysis. The detection of spatial-temporal processes of flooding in satellite data was used to assess the method. A deep-learning approach for change detection applied to satellite time-series images was presented by Sublime and Kalinicheva in [66]. Their approach was compared against other machine-learning methods. The achieved results are better than the results of other existing techniques because of the method’s higher performance and relative fast analysis. An online change detection algorithm was proposed by Chandola and Vatsavai in [78]. They used a Gaussian process based non-parametric time-series prediction model in an online mode. The process solves a large system of equations involving the associated covariance matrix which grows with every time step. Their algorithm identifies changes monitoring the difference between the predictions and previous observations within a statistical control chart framework.

An automatic clearance anomaly detection method was proposed by Chen et al. in [82]. The clearance anomaly detection measures the distance between power lines and other objects (e.g., trees and buildings) to evaluate whether the clearance is within the safe range. The clearance measurements were compared with the standard safe threshold to find the clearance anomalies. The results were validated through qualitatively visual inspection, quantitatively manual measurements in raw point clouds, and on-site field survey. The achieved results show that the proposed method detects the clearance hazards, such as tree encroachment, effectively, and the clearance measurement is accurate.

3.2. Research on Anomaly Detection

Anomaly detection can be used by systems to categorize problems based on the use of taxonomies, such as in [13]. Regarding identification of different types of anomalies useful to help decision-making systems, a unified framework for anomaly detection was proposed by Kittler et al. in [13]. According to them, the scientific literature conventionally understands an anomaly as an outlier. However, their framework has expanded the concept of anomaly beyond the conventional meaning of outlier. They claim that anomalies are related to many other occurrences, such as rare events, unexpected events, distribution drift, noise, and novelty detection of an object or an object primitive. The framework presented the multifaceted nature of anomalies and suggested effective mechanisms to identify and distinguish each of them. The research has provided a taxonomy of anomalies that includes, for example, unknown object, measurement model drift, unknown structure, unexpected structural component, component model drift, and unexpected structure and structural components. Because the practical application of the taxonomy is relatively new, and the detection and identification of the type of the anomaly involves complex aspects, it is appropriate to study the practical application of each type of anomaly individually. The research asserts that by identifying different types of anomalies, systems can select appropriate responses to deal with each type of anomaly. For decision-making systems, the application of the anomaly detection is potentially increased when it is based on incongruence.

3.3. Research on Incongruence

The use of incongruence for detecting the occurrences of some types of events was deeply discussed by Weinshall et al. in [15]. They identified distinct types of events based on the conflicting predictions given by weak and strong classifiers. Their study presented a framework for the representation and processing of those incongruent events. According to them, an incongruent event is an event for which the probability is divergent when computed based on the use of different types of classifiers. They applied their methodology on pictures of dogs and motorcycles, and indoor videos recorded to simulate common situations, such as people walking and talking. Similar to the other studies we have reviewed, this work also does not present the application of their methodology to solve real-world problems.

The decision cognizant Kullback–Leibler (DC-KL) divergence was proposed by Ponti et al. in [16]. The DC-KL divergence reduces the contribution of the minority classes which obscure the true degree of classifier incongruence. Simulations, i.e., synthetic data, were used to analyze the properties of the novel divergence measure. The measure is more robust to minority class clutter, less sensitive to estimation noise, and achieves much better statistics for discriminating between classifier congruence and incongruence than the classical KL divergence.

The delta divergence measure was proposed by Kittler and Zor in [17]. This measure focuses on the most probable hypotheses. Consequently, it reduces the effect of the probability mass distributed over the non-dominant hypotheses. The measure satisfies symmetry and independence of classifier confidence.

A new classifier incongruence measure was proposed by Kittler and Zor in [18]. The measure overcomes the shortcomings of previous measures and presents relatively low sensitivity to estimation noise, under the assumption of constrained Gaussian distribution. Moreover, the measure determines incongruence thresholds at given levels of statistical significance for different measured values corrupted with different levels of noise.

3.4. Research on Real-World Problems, Which Inspired This Study

The studies commented in this subsection inspired us to investigate the potential application of the taxonomy [13] to solve real-world problems. Spatially explicit long-term monitoring frameworks and priority mitigation measures to cope with acute and chronic risks were proposed by Fernandes et al. in [31]. The research focused on the dam failure in Mariana, Brazil, where a dam breach abruptly discharged between 55 and 62 million m³ of ore tailings into the Doce River. They claimed that environmental disasters like that of the Doce River would become more frequent in Brazil. The study is based on differential analyses of Landsat 8 scenes.

An automatic approach to enhance and detect river networks was proposed by Kang Yang in [2]. The research characterizes rivers in accordance with their Gaussian-like cross-sections and longitudinal continuity. River cross-sections were enhanced by a Gabor filter. The approach was applied to detect rivers in Landsat 8 images.

A study which focuses on the mapping and monitoring of the spatial extent of surfaces affected by mine waste was presented by Mielke et al. in [32]. They studied the influence of sensors on the ability to discriminate mine waste surfaces from their surroundings. The research explored the potential use of images from remote sensing to map and monitor the spatial extent of mine waste surface material in areas with mine tailings. They claimed that remote sensing analysis is very important to monitor the extent of mine waste surfaces and that mine waste sites have the potential to contain problematic contaminants, such as chrome, lead, uranium, etc. The research was based on the use of Landsat 8 data.

The impacts of oil extraction on the environment and health of indigenous communities in the Northern Peruvian Amazon (Marañon) were assessed by Rosell-Melé et al. in [33]. They investigated the occurrence of crude oil pollution in soils and river sediments caused by voluntary discharges or accidental spills. Their study is related to a very productive oil area in Peru, which is placed in the headwaters of the Amazon River. Their findings suggest that wildlife and indigenous populations in the Marañon region are exposed to the ingestion of oil, and local spillage of oil in the watercourses could have eventually reached the Amazon River.

4. The Proposed Strategy

4.1. Introduction to the Proposed Strategy

The flowchart in Figure 1 summarizes the methodology of the proposed strategy. More detailed descriptions of this methodology as well as of each parameter used here can be found in Section 5 and Appendix A. All chosen Landsat 8 images (see the fourth column in

Table 1. Information about the studied environments.

Environment and Country	UTM Zone and Center: Latitude/Longitude	Disaster	Landsat 8 Image (Scene)	Date of the: Image Acquisition/Environmental Disaster	Number of Samples
Doce River Brazil	23 20°13′48.07″S 42°43′47.24″W	Dam Breach	LC08_L1TP_217074_20151112_20170402_01_T1	12 November 2015 5 November 2015	250
Athabasca River Canada	12 57°18′35.86″N 112°33′13.21″W	______	LC08_L1TP_043020_20170826_20170913_01_T1	26 August 2017 No Disaster	100
Elbe River Germany	32 54°30′31.36″N 9°29′07.80″E	______	LC08_L1TP_196022_20180509_20180517_01_T1	9 May 2018 No Disaster	100
Marañon River Peru	18 4°20′20.54″S 74°45′45.97″W	Pipeline Disruption	LC08_L1TP_007063_20160827_20170321_01_T1	27 August 2016 22 August 2016	100
Tietê River Brazil	23 23°06′46.33″S 46°30′16.96″W	Permanent Sewage Disposal	LC08_L1TP_219076_20170827_20170914_01_T1	27 August 2017 Permanently	100
Arava Valley Israel	36 30°18′22.86″N 35°02′25.51″E	Pipeline Disruption	LC08_L1TP_174039_20141214_20170416_01_T1	14 December 2014 4 December 2014	100

) were thoroughly evaluated by applying the proposed strategy. Each environment was evaluated separately. Many Landsat 8 scenes were previously checked visually to assess their conditions in order to select only those with low overlay contamination by cloud cover. This is critical to making geographical features as visible as possible. Afterwards, each Landsat 8 scene (see the fourth column in Table 1) was chosen according to its date of acquisition. In cases in which disaster had already occurred, the scene with the date of acquisition closest to the date of the disaster (see the fifth column in Table 1) was chosen in order to best evaluate the environment after the disaster. In the cases in which no environmental disaster had occurred, the most recent images with low overlay contamination by cloud cover were chosen (that is the case of the Athabasca, Elbe and Tietê Rivers).

4.2. The Advantages of Our Strategy

The following are advantages of our strategy: (1) It takes into consideration the multifaceted nature of anomalies exposed by [13]. (2) It helps to offer an appropriate response to deal with anomalies of the type unexpected structure and structural components. (3) To the knowledge of the authors, this is the first effort of a series of studies dedicated to investigate the practical application of detection of anomalies in accordance with the taxonomy published in [13] to solve real-world problems. (4) This research overcomes the number of categories studied in previous investigations, e.g., [25], which presented four different anomalies, because our strategy is based on that taxonomy [13] which presents ten different domain anomalies. (5) It processes just a single Landsat 8 image (scene), whereas other proposed anomaly detection methods use time-series of Landsat 8 images for each environment to be able to detect outliers. Time-series will only be used to detect anomalies of the types measurement model drift and component model drift in future papers.

Additionally, this research has established a relationship between the occurrence of anomalies of the type unexpected structure and structural components and the presence of brown mud in large rivers. Let us consider the example of a mining company which needs to monitor whether the brown mud of its reservoir is polluting a large river. In this example, the mining company does not need to deal with any other type of pollutant other than brown mud, therefore a system able to detect only brown mud is enough for the company. The fewer the images, the fewer the demands for computational resources, consequently making the system cheaper. Moreover, the fewer the images the less data there is to be processed, consequently making the system faster. Therefore, a system based on the application of the proposed strategy, which deals with only one image per environment, will be advantageous. This system will save more money and time for the mining company than systems based on other proposed anomaly detection methods, because they generally deal with time-series.

For the case studies presented in this paper, the anomalies are related to brown mud or crude oil present as either potential contaminators or contaminators of large rivers. Our strategy identifies incongruences to detect oil spills in some areas where there is water, or brown mud, and where there are high levels of turbidity in the water. These findings are important because an incongruence-based anomaly detection strategy can greatly increase the ability of surveillance systems to detect environmental disasters and to perform mappings.

5. Materials and Method

5.1. Materials

Landsat 8 images (scenes of environments) from the United States Geological Survey (USGS) data sets and benchmarks [83] were used to evaluate the proposed strategy. The USGS consists of Landsat images and many other important data set collections, which are all freely available in [83] for non-commercial purposes. These data sets were chosen to perform the experiments in this study because they provide diversity regarding the images and geographical features for many different environments. Moreover, they have been frequently used to evaluate many recently-published proposals, such as [2,31,32], which allows us to compare the present proposed strategy with other studies. For this study, each scene was provided by the USGS through the Earth Explorer Platform.

All images used in this study are qualified with value “9,” which represents the best quality for Landsat 8 images [26]. Moreover, all Landsat 8 data used in the present experiments were collected by the operational land imager (OLI), which is an instrument onboard the Landsat 8 satellite. The size of each image is approximately 170 km north-south by 183 km east-west, i.e., 106 mi by 114 mi. These images are available in the tagged image file format (TIFF). The three most common types of TIFF images are gray scale, indexed, and RGB.

For this study, each task of our strategy, including all classifications, was systematically performed using the Qgis 2.18.19 (Las Palmas) software with the Orfeo tool box. The software and the tool box are all freely available in [84] for non-commercial purposes. Qgis is a well-known and extensively used software composed of many tool sets, which are useful for analyzing images, mainly from remote sensing, for geosciences.

5.2. Study Areas

Figure 2 shows the scenes and the hydrographical maps of the six different areas studied in this research. The information about the studied environments, including the identification of the Landsat 8 images (scenes) used in the experiments, is listed in Table 1. The geographical settings are listed in

Table 2. Geographical settings.

Environment and Country	Basin Size	Length of the River	Predominant Features in the Environment	Human Activity in the Region	Potential Threat to the Environment	Potential Contaminators or Contaminators
Doce River Brazil	83,400 km2 (32,201 sq mi)	853 km (530 mi)	Mountains and Valleys	Mining	Ore Tailings Reservoirs	Brown Mud
Athabasca River Canada	95,300 km2 (36,800 sq mi)	1231 km (765 mi)	Lakes	Mining	Ore Tailings Reservoirs	Brown and Red Mud
Elbe River Germany	148,268 km2 (57,247 sq mi)	1094 km (680 mi)	Close to the Coastline	Mining	Ore Tailings Reservoirs	Red Mud
Marañon River Peru	358,000 km2 (138,000 sq mi)	1737 km (1079 mi)	Forest	Transport of Crude Oil	Oil Pipeline	Crude Oil
Tietê River Brazil	150,000 km2 (58,000 sq mi)	1150 km (710 mi)	Densely Populated Region	Sewage Disposal	Sewer	Sewage
Arava Valley Israel	______	______	Arid Valley	Transport of Crude Oil	Oil Pipeline	Crude Oil

. All data are related to areas of large river basins, except the area of the Arava Valley. The image of the Arava Valley in southern Israel, acquired by the Landsat 8 satellite on 14 December 2014, was chosen because its soil received an estimated 3–5 million liters of crude oil on 4 December 2014 (29°40′29.28″N 35°0′25.56″E) [85]. This image was used to contrast the experiments with the Marañon River since there is absence of water in that part of the Arava Valley. This contrast is important because only two classes were defined for all experiments: “water,” which categorizes any open water body, and “no-water,” which categorizes anything different from open water bodies, for example land, trees, cities (this approach is similar to the one found in [86]). The idea is to evaluate if the proposed strategy would detect oil as if it were water. If oil were to be detected as water, then all experiments with images of the Marañon River basin would be invalidated. The image of the Marañon River in Peru, acquired by the Landsat 8 satellite on 27 August 2016, was chosen because its water received high quantity of crude oil on 22 August 2016. This oil spill was the consequence of a Colombian environmental disaster caused by a pipeline disruption (4°48′45.6″S 75°23′56.8″W) [87]. The crude oil is present around almost all extensions of the same pipeline since it was the source of other similar environmental disasters in the preceding three years.

The image of the Doce River in Brazil, acquired by the Landsat 8 satellite on 12 November 2015, was chosen because its water received around 55–62 million m³ of iron ore tailings (brown mud) on 5 November 2015, as a result of Brazil’s worst environmental disaster caused by a dam breach (20°12′23.4″S 43°28′01.6″W) [31]. The image of the Elbe River in Germany, acquired by the Landsat 8 satellite on 9 May 2018, was chosen because of the location of a dam close to its water (53°38′20.4″N 9°25′02.0″E). The river could receive many millions of m³ of bauxite tailings (red mud) in the case of an environmental disaster similar to that which occurred in Brazil. The image of the Athabasca River in Canada, acquired by the Landsat 8 satellite on 26 August 2017, was chosen because there is a set of dams located close to its water (57°03′06.3″N 111°35′18.0″W), which could cause the river to receive many millions of m³ of brown or red mud in the case of an environmental disaster similar to that which occurred in Brazil. The image of the Tietê River in Brazil, acquired by the Landsat 8 satellite on 27 August 2017, was chosen because its water receives significant quantities of sewage each year (23°23′24.3″S 46°58′36.5″W).

For each Landsat 8 image used in this study, at least 100 samples were located and validated based on the topographic maps of the environments, adopted as ground truth for the purpose of this work. There are samples spread all over the scene. The numbers of samples used as training data to model each environment according to this approach are shown in the sixth column of Table 1, being 50% constituted by open water bodies and 50% by other geographical features. Regarding samples of open water bodies which are expected to be classified as “water” by the classifiers, the chosen areas were carefully selected from lakes, small rivers, artificial reservoirs, waterways, ocean (if present), etc., present in the image. Regarding samples of other geographical features which are expected to be classified as “no-water” by the classifiers, the chosen areas were carefully selected from fields, forests, clouds, cities, highways, mountains, plantations, shadows of clouds, etc., present in the image.

Each image presents different difficulty levels to locate and analyze geographical features, in the environments. For example, in any of these environments it is possible to find parts of the image with an easy to detect, large, and singular geographical feature that differs from other parts of the image with hard to detect, small, numerous, and geographically spread features. By using diversified environments, it was possible to verify the accuracy of the proposed strategy to identify incongruences in order to detect anomalies, when it is applied to different scenes as well as to different types of geographical features present in images.

5.3. Method

5.3.1. Data Preprocessing

Each image was processed performing tasks to identify incongruences in accordance with Figure 1. Initially, seven bands from a Landsat 8 scene were added as raster layers to a new project of Qgis: Band 1—coastal aerosol (0.43–0.45 µm) 30 m; Band 2—blue (0.45–0.51 µm) 30 m; Band 3—green (0.53–0.59 µm) 30 m; Band 4—red (0.64–0.67 µm) 30 m; Band 5—near infrared (NIR) (0.85–0.88 µm) 30 m; Band 6—SWIR 1 (1.57–1.65 µm) 30 m; and Band 7—SWIR 2 (2.11–2.29 µm) 30 m. This first task was done in order to organize the seven aforementioned bands into a single multiband raster to allow easy access.

Second, a band composition was built in order to allow the user of the system to analyze the scene visually. A band composition is an image composed of different bands of a satellite image. The band composition uses pixels of different color values in order to represent various properties of objects from the real world. A band composition is represented by Equation (1), for which, considering each Landsat 8 band as a two-dimensional array of pixels, x and y are coordinates of the pixels. The composition vr is a tridimensional array composed of three overlapped elements (bands), being b4(x, y) = Band 4, b5(x, y) = Band 5, and b6(x, y) = Band 6. Figure 3 shows the band composition R(4)G(5)B(6) as an example.

$$vr\left(x,y\right)=\left[\begin{array}{c}b4\left(x,y\right)\\ b5\left(x,y\right)\\ b6\left(x,y\right)\end{array}\right]$$

(1)

Next, we applied a histogram stretching to the composite R(4)G(5)B(6). The procedure was performed by a Qgis tool from the menu properties (style). Because Qgis needs to create a color table in order to render the band composite, the tool calculates the table based on the mean and standard deviation for the three selected bands.

The mean and the standard deviation were calculated respectively according to Equations (2) [88] and (3) [88], for which, x and y are coordinates of the pixels, g_i(x, y) is an image, K is the total number of images g(x, y), i is the identifier of each image g(x, y), and f(x, y) is the image formed averaging K images g(x, y). In Equation (3), h(x, y) is the image formed calculating the standard deviation. Figure 4 shows an example of an enhanced image [89].

$$f\left(x,y\right)=\frac{1}{K}{\displaystyle \sum }_{i=1}^K{g}_i\left(x,y\right)$$

(2)

$$h\left(x,y\right)=\sqrt{\frac{1}{K}{\displaystyle \sum }_{i=1}^K {\left(g_i\left(x,y\right)-f\left(x,y\right)\right)}^2}$$

(3)

Shortly after, we created a new project in Qgis to receive the Landsat 8 processed image. The chosen coordinate reference system (CRS) was WGS84 and with the universal transverse mercator (UTM) map projection (see the second column in Table 1). Another band from the Landsat 8 scene was added as a raster layer to the project of Qgis: Band 8 panchromatic (PAN) (0.50–0.68 µm) 15 m. Next, a pan-sharpening step based on component substitution (CS) or projection substitution was performed [89]. Equation (4) represents a general formulation of this CS method [90],

$${\widehat{MS}}_k={\widetilde{MS}}_k+g_k\left(P-I_L\right),\hspace{1em}\hspace{1em}k=1,2,\dots ,N$$

(4)

in which MS is a multispectral image, $\widehat{MS}$ is a pan-sharpened image, $\widetilde{MS}$ is a multispectral image interpolated at the scale of the panchromatic image, the subscript k indicates the k^th spectral band, g = [g₁,…, g_k,….,g_N] is the vector of the injection gains, P is the histogram-matched panchromatic image, and I_L is defined by the Equation (5) [91],

$$I_L={\displaystyle \sum }_{i=1}^N{w}_i{\widetilde{MS}}_i$$

(5)

in which the weight vector w = [w₁,…, w_i,…, w_N] is the first row of the forward transformation matrix and may be chosen, whenever possible, to measure the degrees of spectral overlap among the multispectral and panchromatic channels.

Figure 5 presents a flowchart summarizing how the pan-sharpening method based on this CS method works. The flowchart has blocks responsible to interpolate the multi-scale image for matching the scale of the panchromatic image, calculate the intensity component based on Equation (5), match the histograms of the panchromatic image and the intensity component, and inject the extracted details based on Equation (4).

For the pan-sharpening method, interpolation must guarantee the overlap of multispectral and panchromatic images at a finer scale. Depending on the acquisition geometry of the imaging instruments, interpolation with conventional zero-phase linear finite-impulse response filters may require a manual realignment, e.g., after a bi-cubic interpolation. Alternatively, linear nonzero-phase filters, having even numbers of coefficients, may be profitably used, according to [91,92].

Although the 30-meter original spectral bands are enough for dealing with large rivers, the pan-sharpening was necessary for dealing with lanes that support oil pipelines. Those lanes, although long, are much narrower than large rivers. Working with 15-meter images made possible the detection of the presence of oil on those narrow lanes. It is important because, for example, there is a very long lane, built in the vicinity of the Marañon River in order to support oil pipelines, that has already received high quantities of crude oil as a consequence of successive oil spills. Whereas those lanes should normally present either dry or humid soil around the pipelines, a significant oil spill turns the soil impermeable causing accumulation of water around the pipelines. The detection of water retained in the lane, therefore, is indicative of the occurrence of an oil spill on the soil of that area.

Next, the pan-sharpening resulting image had its histogram stretched. This step was performed exactly the way it was previously mentioned, including the calculation of Equations (2) [88] and (3) [88].

Then, the second-order image statistics were calculated using the Orfeo toolbox in order to compute the global mean and standard deviation for each band from a set of images. Their applications are related to geometric modeling, i.e., to create a model taking into consideration the spatial distribution of the pixels in the image. The model was created in order to statistically represent the enhanced image, and the statistics are not affected by the enhancing parameters. The second-order image statistics are those for which the slope of the power spectrum tends to be close to negative two. The power spectrum of an M by M image is represented by Equation (6) [93], for which, F is the Fourier transform of the image, φ are directions, and u and v are two-dimensional frequencies represented in polar coordinates in accordance with Equations (7) and (8) respectively.

$$S\left(u,v\right)=\frac{{\left|F\left(u,v\right)\right|}^2}{M^2}$$

(6)

$$u=f\cos \varnothing$$

(7)

$$v=f\sin \varnothing$$

(8)

According to Reinhard et al. in [93]:

“(...) averaging over all directions φ and all images in the ensemble, it is found that on log-log scale power as function of frequency f lies approximately on a straight line (...) This means that spectral power as function of spatial frequency behaves according to a power law function. Moreover, fitting a line through the data points yields a slope α of approximately negative two for natural images:

$$S\left(f\right)\propto A{f}^2=A{f}^{-2-n}$$

(9)

Here, α ≈ −2 is the spectral slope, η is its deviation from −2, and constant A describes the overall image contrast.”

Contrast can be obtained by image statistics, for example, calculating the standard deviation of all pixel intensities divided by the mean intensity (σ/µ) in accordance with Equation (10).

$$\frac{{\sigma }^2}{{\mu }^2}={\displaystyle \sum }_{\left(u,v\right)}S\left(u,v\right)$$

(10)

Afterward, in the image, at least 100 samples were selected and validated based on maps of the environments (ground truth). The details regarding the chosen samples were described in the Section 5.2. The selection of samples is important because they are used as training data to model each environment for contextual or non-contextual classifiers. Whenever users select samples for classifiers, their classifications are named supervised or semi-automatic. A supervised classification allows the user to collect samples from the image and to use a classifier trained with these samples to classify the image [90].

The preprocessing is essential to prepare the data to be used for the learning and classification approaches adequately. Our main contribution in the preprocessing is regarding the care for selecting the samples. The samples were not located on the large rivers, ore tailings reservoirs, oil pipelines, or lanes for oil pipelines, in order to avoid the experiments yielding biased results. In other words, no areas located on the large rivers were selected to be samples because these samples could bias the classifiers to classify the large rivers as “water.” No areas located on ore tailings reservoirs, oil pipelines, or lanes for oil pipelines were selected to be samples because these samples could bias the classifiers to classify those structures as “no-water.”

5.3.2. Learning and Classification Approaches

The training of classifiers step was performed based on the enhanced 15-meters image, the statistical model, and the selected samples. The training is important in order to create the models which will be used by the contextual or non-contextual classifiers. The classifiers use these models as references based on the selected samples to classify features of interest present in the enhanced 15-meters image [94]. Next, classifications were applied using various classifiers, such as Naïve Bayes [35], super vector machine (SVM) [8,65], decision tree (DT) [29,30], k-nearest-neighbor (kNN) [27,28], and boost [9]. Among the tested classifiers, the best results were achieved by kNN and decision tree (non-contextual classifiers), and boost (contextual classifier). Although the other tested classifiers were used the same way as kNN, DT, and boost, they were excessively time-consuming. Moreover, their results did not overcome those achieved by kNN, DT, and boost for the experiments of this study. Therefore, we decided to consider only kNN, DT, and boost as the classifiers of the experiments in this study.

The kNN classifier is useful to perform discriminant analysis in situations for which it is hard to determine parametric estimates of probabilities [27]. Let us consider that X = [x₁,...,x_N] is the training data with N points of dimensionality D, X_i = [x_i1,...,x_ik] is the k nearest neighbors of x_i, X_t is the testing data with N_t points, x₀ is an arbitrary testing data point, X₀ = [x₀₁,...,x_0k] contains its k nearest neighbors from training data, for which the labels are [l₁,...l_k], and Ω = [Ω₁,…, Ω_C] is the set composed of C classes present in the data. Additionally, the kNN considers that the k neighbors of a testing point have equal weights [28].

The classification performed based on the kNN finds the nearest neighbors regarding a testing point in the training data. Moreover, the classifier assigns the test point to the most frequently occurring class of its k neighbors. The kNN classifies x₀ based on the majority voting rule presented by Equation (11) [28], for which δ is the Kronecker delta represented by Equation (12).

$$j^{*}=arg\begin{array}{c}max\\ j=1,\dots ,C\end{array} \left({\displaystyle \sum }_{i=1}^k\delta \left(l_i,j\right)\right)$$

(11)

$${\delta }_{ij}=\{\begin{array}{c}1,\hspace{1em}ifi=j\\ 0,\hspace{1em}ifi\ne j\end{array}$$

(12)

The decision tree classifier is useful to perform decision analysis by finding the most probable decision for achieving an objective [29]. Many algorithms are able to build decision trees. Among them, the Iterative Dichotomiser 3 (ID3) is a well-known algorithm used to generate decision trees from data sets. The ID3 is based on information entropy and information gain, which are represented respectively by Equations (13) and (14) [30].

$$Ent\left(D\right)=-{\displaystyle \sum }_{d=1}^k{p}_k{\log }_2{p}_k$$

(13)

$$Gain\left(D,a\right)=Ent\left(D\right)-{\displaystyle \sum }_{v=1}^V\frac{\left|D_v\right|}{\left|D\right|}Ent\left(D_v\right)$$

(14)

These equations consider that the training sample set is represented by D = {(x₁, y₁), (x₂, y₂),…, (x_m, y_m)}, and |D| stands for the training sample number. The p_k (k = 1, 2, …, |D|) is the ratio of every type for every current sample set. A = {a₁, a₂ · · · a_d} stands for the attribute set of D, e.g., color, shape, texture, etc., for which d = {1, 2, …, k}. For each attribute a_i, V is a set of features composed of v different values and represented by V = {$a_i^1$, $a_i^2$, …, $a_i^v$}, e.g., $a_1^1$ = green, $a_1^2$ = blue, etc. D_v is the sample subset related to the value $a_i^v$ from a_i in D, and |D_v| stands for the number of the current samples of the subset. The entropy of D, Ent(D), characterizes the impurity of the collection of samples in D, i.e., it measures the amount of uncertainty in the data set D. The gain of D, Gain(D,a), is the expected reduction in the entropy of D related to the choice of the attribute a.

The Boost classifier combines iteratively weak classifiers by taking into consideration a weight distribution on the training samples such that more weight is attributed to samples misclassified by the previous iterations [9]. The final strong classifier is a weighted combination of weak classifiers followed by a threshold. The steps performed by the Boost are described in detail as follows.

(1) Establish training sets.

Let us consider that x_i and y_i are respectively the samples and their labels. Moreover, D_t(x_i,y_i) is the t^th distribution over all these training samples. Let us also consider a training set {(x₁,y₁),(x₂,y₂),...,(x_n,y_n)}, for which y_i ∈ {-1,1}.

(2) Initialize weights.

The initialization of weights is performed according to Equation (15) [9].

$$D_0\left(x_i,y_i\right)=\frac{1}{N},\hspace{1em}\hspace{1em}i=1,2,\dots ,N$$

(15)

(3) Training loop.

For t = 1,2,...,T (T is the maximum training number):

(a) For each feature j, train a classifier h_j which is a simple linear classifier, i.e., a classifier restricted to use a single feature. Equation (16) expresses this type of classifier [9].

$$h\left(x_i\right)=\{\begin{array}{ll}1,& if p_j{x}_{i,j}<p_j{\theta }_{i,j}\\ -1,& otherwise\end{array}$$

(16)

In Equation (16), x_i,j is the value of the j^th feature of the sample x_i, p_j ∈ {−1,1} determines the direction of the inequality sign, θ_i,j denotes the threshold value of the j^th feature of the sample x_i. Equation (17) shows that the error ε_t is evaluated with respect to D_t(x_i,y_i) [9].

$${\epsilon }_j={\displaystyle \sum }_{i:y_i\ne h_j}{D}_t\left(x_i,y_i\right)$$

(17)

(b) Choose the weak classifier h_t, with the lowest error ε_t.

(c) If ε_t ≥ 1/2, then stop.

(d) If ε_t < 1/2, calculate α_t which is the weight assigned to the classifier h_t, in accordance with Equation (18).

$${\alpha }_t=\frac{1}{2}In\left(\frac{1-{\epsilon }_t}{{\epsilon }_t}\right)$$

(18)

(e) Equation (19) shows how to update the weights [9],

$$D_{t+1}\left(x_i,y_i\right)= \frac{D_t\left(x_i,y_i\right)e^{-{\alpha }_t{y}_t{h}_t\left(x_i\right)}}{Z_t}$$

(19)

for which Z_t is a normalized constant computed to insure that D_t(x_i,y_i) represents a true distribution, established by Equation (20) [9].

$${\displaystyle \sum }_{i=1}^N{D}_t\left(x_i,y\right)=1$$

(20)

(4) Output the final classifier in accordance with Equation (21) [9].

$$H\left(x\right)=sign \left({{\displaystyle \sum }}_{t=1}^{Tk}{\alpha }_t{h}_t\left(x\right)\right)\hspace{1em}\hspace{1em}\hspace{1em}sign\left(x\right)=\{\begin{array}{l}-1, x<0\\ 0, x=0\\ 1, x>0\end{array}$$

(21)

The image classifications are important in order to group features of interest into classes according to the similarities of their characteristics. In the images resulting from the classifications, each class is represented by a different color. In the case of this research, because only two different classes were used for grouping features of interest, “water” and “no-water,” the images resulting from the classification processes are binary [94].

At the end, a subtraction of images resultant from non-contextual and contextual classifications was performed, computing the difference of all pairs of corresponding pixels from both images [88]. The subtraction was calculated according to Equation (22) [88], for which, x and y are coordinates of the pixels, and f(x, y) and h(x, y) are respectively the images resultant from non-contextual and contextual classifications. Image g(x, y) results from the subtraction. This subtraction is important because it allows the automatic identification of congruence, in which case the subtraction results in zero, or automatic identification of incongruence, in which it results in at least one pixel equal to one. Incongruences can reveal the presence of anomalies as explained in [16]. When anomalies were identified, they were analyzed and categorized by type in accordance with [13].

$$g\left(x,y\right)=f\left(x,y\right)-h\left(x,y\right)$$

(22)

The learning and classification approaches are essential to make feasible the application of the taxonomy [13] to solve real-world problems. Our main contributions in the learning and classification approaches are as follows: (1) The use of a single Landsat 8 image previously assessed as high-quality sensory data. (2) The creation of statistical and classification models adequate enough to make PR tools learn how to classify the features present in the environments. (3) The subtraction of the results from contextual and non-contextual classifiers to indicate the occurrence of incongruence. (4) Finding the combination of the three previous conditions to establish a relation between the occurrence of anomalies of the type unexpected structure and structural components and the presence of brown mud in large rivers.

Regarding methodological limitations, it is recognized that the higher the resolution of the images, the longer the computational time to process the images. Therefore, long time of processing is inevitable when large images are evaluated by anomaly detection strategy. This time depends on the computer being used and nowadays can vary from dozens of minutes to many hours. However, this variation, in the computational time, spent by the proposed strategy to process the images is acceptable since different computers have different computational resources.

6. Results

Three evaluation measure mechanisms (contextual and non-contextual classifiers; sensory data quality assessment; incongruence indicator) were used to reliably qualify the anomaly, as described in [13], because they are three successful mechanisms to qualify anomalies. According to Kittler et al. [13], three conditions need to be satisfied to guarantee the occurrence of the anomaly of the type unexpected structure and structural components. The three conditions are high sensory data quality, contextual and non-contextual classification of the same image, and incongruence. If any of the three conditions are not present, it is not possible to identify an anomaly of this type. In this study, the high sensory data quality was confirmed because all images are qualified with value “9,” which represents the best quality for Landsat 8 images. The contextual and non-contextual classifications were confirmed by the results achieved applying respectively boost and kNN or DT classifiers. The presence of incongruence results was confirmed performing the last step of the strategy.

All assessments were carefully performed taking into consideration the application of the proposed strategy on Landsat 8 images. Since all those Landsat 8 images are large (170 km north-south by 183 km east-west, i.e., 106 mi by 114 mi, approximately), the image of each environment was precisely cropped in order to create a set of smaller images tiles. This was fitting, as the majority of assessment approaches found in the scientific literature are applied to smaller images.

Table 3. Information about how cropping was done for each Landsat 8 scene.

Environment Country	Landsat 8 Scene	Scene Resolution (Height and Width in Pixels)	Image Tile Resolution (Height and Width in Pixels)	Total Number of Tiles
Doce River Brazil	LC08_L1TP_217074_20151112_20170402_01_T1	15705 × 15440	151 × 193	8400
Athabasca River Canada	LC08_L1TP_043020_20170826_20170913_01_T1	16499 × 16277	146 × 155	11836
Elbe River Germany	LC08_L1TP_196022_20180509_20180517_01_T1	16099 × 15882	145 × 143	12544
Marañon River Peru	LC08_L1TP_007063_20160827_20170321_01_T1	15498 × 15185	126 × 146	12792
Tietê River Brazil	LC08_L1TP_219076_20170827_20170914_01_T1	15484 × 15250	139 × 162	14641
Arava Valley Israel	LC08_L1TP_174039_20141214_20170416_01_T1	15432 × 15067	129 × 127	10528

shows how the cropping was done for each Landsat 8 scene.

In order to quantitatively validate the proposed strategy, we used the metrics accuracy, precision, recall [15], and F-measure [82] for which, TP are true positives, i.e., images in which incongruences are truly detected; TN are true negatives, i.e., images in which congruences are truly detected; FN are false negatives, i.e., images in which congruences are falsely detected; FP are false positives, i.e., images in which incongruences are falsely detected. Accuracy measures the efficiency of results and is represented by Equation (23), for which M is the number of images. Precision measures the relevancy of results and is represented by Equation (24). Recall measures the quantity of truly relevant results and is represented by Equation (25). F-Measure measures the balance between precision and recall and is represented by Equation (26).

$$Accuracy=\frac{TP+TN}{M}$$

(23)

$$Precision=\frac{TP}{TP+FP}$$

(24)

$$Recall=\frac{TP}{TP+FN}$$

(25)

$$F-Measure=2\times \frac{Precision\times Recall}{Precision+Recall}$$

(26)

6.1. Experimental Results

Table 4. Contingency table helping to detect incongruences in the Landsat 8 scene of Doce River consisting of 8400 tiles (TP—true positive, FP—false positive, FN—false negative, TN—true negative).

	Incongruent Event	Congruent Event
Detection Incongruent	TP = 68	FP = 20
Detection Congruent	FN = 0	TN = 8312

,

Table 5. Contingency table helping to detect incongruences in the Landsat 8 scene of Athabasca River consisting of 12,084 tiles (TP—true positive, FP—false positive, FN—false negative, TN—true negative).

	Incongruent Event	Congruent Event
Detection Incongruent	TP = 87	FP = 38
Detection Congruent	FN = 0	TN = 11,959

,

Table 6. Contingency table helping to detect incongruences in the Landsat 8 scene of Elbe River consisting of 12,544 tiles (TP—true positive, FP—false positive, FN—false negative, TN—true negative).

	Incongruent Event	Congruent Event
Detection Incongruent	TP = 0	FP = 18
Detection Congruent	FN = 0	TN = 12,526

,

Table 7. Contingency table helping to detect incongruences in the Landsat 8 scene of Marañon River consisting of 12,915 tiles (TP—true positive, FP—false positive, FN—false negative, TN—true negative).

	Incongruent Event	Congruent Event
Detection Incongruent	TP = 39	FP = 13
Detection Congruent	FN = 0	TN = 12,863

,

Table 8. Contingency table helping to detect incongruences in the Landsat 8 scene of Arava Valley consisting of 10,528 tiles (TP—true positive, FP—false positive, FN—false negative, TN—true negative).

	Incongruent Event	Congruent Event
Detection Incongruent	TP = 0	FP = 25
Detection Congruent	FN = 0	TN = 10,503

and

Table 9. Contingency table helping to detect incongruences in the Landsat 8 scene of Tietê River consisting of 14,641 tiles (TP—true positive, FP—false positive, FN—false negative, TN—true negative).

	Incongruent Event	Congruent Event
Detection Incongruent	TP = 0	FP = 0
Detection Congruent	FN = 0	TN = 14,641

, which are all contingency tables [95,96,97], present the quantitative evaluations of incongruence detections on real images for, respectively: Doce River, Athabasca River, Elbe River, Marañon River, Arava Valley, and Tietê River.

6.2. Interpretation of the Results

Table 4, Table 5, Table 6, Table 7, Table 8 and Table 9 present encouraging results. As can be seen, there is a significant quantity of true positives associated with the image sets that present brown mud, whereas this quantity is slightly lower for images presenting crude oil and considerably lower for images presenting red mud and sewage. False positives are present in almost all image sets, but their number are of little significance for this study.

The results regarding the accuracy, precision, recall, and F-measure of our strategy are shown in

Table 10. Accuracy, precision, recall, and F-measure of the proposed strategy for different environments (Landsat 8 scenes).

Landsat 8 Scene	Environment	Accuracy	Precision	Recall	F-Measure
LC08_L1TP_217074_20151112_20170402_01_T1	Doce River Basin	~99.76%	~77.27%	100%	87.18%
LC08_L1TP_043020_20170826_20170913_01_T1	Athabasca River Basin	~99.69%	~69.6%	100%	82.08%
LC08_L1TP_007063_20160827_20170321_01_T1	Marañon River Basin	~99.90%	75%	100%	85.71%

.

Table 11. Comparison of the accuracy, precision, recall, and F-measure of this study with others.

Study	Accuracy	Precision	Recall	F-Measure
Ours	~99.78%	~73.96%	100.00%	85.04%
[98]	99.59%	----------	----------	----------
[81]	99.20%	91.85%	53.55%	67.66%
[99]	99.14%	----------	----------	----------
[76]	98.49%	83.84%	83.66%	83.76%
[80]	98.00%	----------	----------	----------
[25]	~91.20%	98.10%	95.7%	96.88%
[77]	88.68%	90.62%	79.62%	84.76%
[66]	84.00%	63.00%	81.00%	70.88%
[100]	81.18%	----------	----------	----------
[78]	78.00%	82.00%	75.00%	78.34%
[82]	----------	96.50%	94.8%	95.64%

shows the comparison of the accuracy, precision, recall, and F-measure of this study, presented in the first row, with results achieved by other studies. It seems that the results are quantitatively consistent with those found in the scientific literature. On average, our strategy achieved the highest values of accuracy and recall among the presented studies. Although our strategy did not achieve the highest values of precision and F-measure compared to the other studies, it is of little significance because the achieved precision and F-measure are high.

Our strategy achieved the highest value of F-measure among the presented studies, except for [25] and [82]. However, [25] and [82] are very specific case studies. Whereas in this paper and in [66,76,77,78,80,81,82,98,99,100] the areas of study are heterogeneous, in [25], results were presented only related to an ice shelf in the Antarctic region, which is a homogeneous area of study. According to [100], results tend to be higher for homogeneous areas of study compared to heterogeneous ones. Whereas in this paper and in [66,76,77,78,80,81,82,98,99] the altitude of the remote sensing is high; in [82], the sensor was very close to the investigated outliers when the images were captured. The unmanned aerial vehicle (UAV) which performed the remote sensing flew in a short distance over the area of study. The higher the proximity of the area of study, the higher tend to be the precision and recall of the detection, consequently making the F-measure higher. Therefore, our strategy achieved the highest value of F-measure among the presented studies, for which the characteristic of the studied area is heterogeneous and the altitude of the remote sensing is high.

Figure 6, Figure 7 and Figure 8 show the examples of the application of the proposed strategy on three different images (tiles) taken respectively from the Doce River, the Athabasca River, and the Marañon River scenes: (a) enhanced image, (b) result overlay on the real image, (c) contextual classifier output, and (d) non-contextual classifier output. The incongruences can be recognized comparing the differences between the outputs (c) and (d) of each image.

In the first case, apart from other geographical features, Figure 6a shows the non-polluted water of a narrow and long lake, in black, in the bottom right corner and the water of part of the Doce River polluted with ore tailing waste (brown mud), also in black, starting on the left, making a sharp turn, and ending on the top. Figure 6c shows the capability of the contextual classifier (Boost) for recognizing both polluted and non-polluted water, whereas Figure 6d shows the capability of the non-contextual classifier (decision tree) for recognizing only non-polluted water.

In the second case, apart from other geographical features, Figure 7a shows part of a large reservoir polluted with ore tailing waste (brown mud), in black, in the top right corner and the non-polluted water of part of the Athabasca River, also in black, starting on the bottom left, making a sharp turn, and ending on the top. Figure 7c shows the capability of the contextual classifier (boost) for recognizing both polluted and non-polluted water, whereas Figure 7d shows the capability of the non-contextual classifier (decision tree) for recognizing only non-polluted water.

In the third case, apart from other geographical features, Figure 8a shows part of a very long lane built in the vicinity of the Marañon River in order to support oil pipelines, in black and orange, starting on the top left and ending on the top right. Figure 8c shows the capability of the contextual classifier (Boost) for recognizing the polluted water concentrated on the impermeable soil in the lane, whereas Figure 8d shows that the non-contextual classifier (kNN) is not so much capable of recognizing the polluted water in the lane.

In the last step of our method, the image of the Figure 6d, Figure 7d and Figure 8d were subtracted respectively from the images of the Figure 6c, Figure 7c and Figure 8c. Figure 6b, Figure 7b and Figure 8b show the respective results of these subtractions highlighted in white and overlapped on the real images. The small white points spread on Figure 6b and Figure 7b can be easily removed applying a morphological opening filter in the image. However, the same filter would remove the result present in Figure 8b. Therefore, we have decided to not apply the filter.

7. Discussion

When the classes are very imbalanced, precision and recall measure the success of prediction. If the classifiers return correct results and the majority of all positive results, then precision and recall present high values. We achieved high values of accuracy, precision, and recall with the application of the incongruence based anomaly detection strategy in this study. Therefore, a significant quantity of truly relevant results showed that the anomalies of the type unexpected structure and structural components were detected with efficiency by our strategy. The values of precision would be even higher if we had limited the studied areas to the vicinity of the locations with pollutants instead of studying the whole Landsat 8 scene, since many false positives were located far from these locations. Moreover, the achieved results demonstrate that the application of the proposed strategy to detect water pollution is consistent with other studies, although this study has introduced an approach to reach results that is unlike the approaches commonly used.

The results presented in the first row of Table 10 are aligned with the results reached by [31], which reveal high concentration of brown mud in the Doce River. Moreover, the results presented in the second row of Table 10 compare well with the results reached by [32], which highlight the ability to discriminate mineral mine waste surfaces. Therefore, the present achieved results suggest that the incongruence occurs because, contrary to the non-contextual classifier (decision tree), the contextual one (boost) is able to detect the water, even for cases in which there is high turbidity in the water of the studied river or reservoirs, because of the high concentration of brown mud.

The results presented in the third row of Table 10 are parallel to the results reached by [87], which also correspond to detection of oil spills in optical satellite images. In the case of our study, instead of detecting oil spills on the surface of the water such as in [87], the detection was performed in lanes built in the vicinity of the Marañon River in order to support oil pipelines. A significant oil spill turns the soil impermeable causing accumulation of water around the pipelines, whereas those lanes should normally present either dry or just humid soil around the pipelines. Therefore, the present achieved results suggest that the incongruence occurs because, contrary to the non-contextual classifier (kNN), the contextual one (Boost) presents high sensibility to detect the water around the pipelines (see Figure 9), revealing the presence of oil spill on the soil of that area. Additionally, the results presented in Table 8, which are related to the Arava Valley, have demonstrated that neither the non-contextual nor the contextual classifier classify crude oil as water in a dry environment, which corroborates the causes of the incongruence previously mentioned.

Outcomes from this study align with those presented by Kittler et al. in [13] and Weinshall et al. in [15], because they perform research involving the use of incongruences and recognize their potential to be applied to anomaly detection. Additionally, the current findings expand these prior studies to real-world images from remote sensing, instead of only synthetic data, and apply to real-world applications, such as the detection of water pollution. For example, after detecting anomalies, decision-making systems can guide computers to give adequate and fast responses to prevent water pollution as an environmental disaster or to analyze the extent to which they occur. Furthermore, the proposed strategy broadens a list of studies, such as [55,56,57,58,59,60,101], engaged to analyze images from remote sensing to offer possible solutions to detect and monitor water pollution.

The results achieved by this study demonstrate higher accuracy and recall than previous studies’ results, such as [25,66,76,77,78,80,81,82,98,99,100]. Therefore, the results suggest that our strategy provides more efficient results and a higher quantity of truly relevant results than the referred studies. The smallest difference regarding accuracy in comparison to our study was reported by Che et al. in [98]. The high level of accuracy was achieved by applying a method based on machine learning. However, Che et al. did not report values in regard to recall.

Another small difference regarding accuracy in comparison to our study was reported by Bormann et al. in [81]. In that study, a method based on detection index allowed them to achieve a high level of accuracy. In contrast, [81] reported the largest difference regarding recall in comparison to our study. Therefore, our method provides much more truly relevant results than the one proposed by them.

Three other studies [76,80,99] reported accuracies, whose values are lower than the one achieved by our study, although they are still close. In [76], a method based on median absolute deviation provided high level of accuracy. However, our method provides more truly relevant results than the one proposed in [76]. Similarly to our study, SVM was a method applied for classification by Bhaduri et al. in [80]. In contrast, we achieved the best results by applying other classifiers in our study compared to the use of SVM. A significant result was provided by a deep-learning method based on convolutional neural network (CNN) in [99]. Although [80,99] achieved high levels of accuracy applying their methods, they did not report results in regard to recall.

Other studies, such as [66,77,78,100], have more significant differences regarding both accuracy and recall in comparison to our study. In these studies, the levels of accuracy and recall were achieved after applying a method based on: (1) seasonal autoregressive integrated moving average (SARIMA) model for autocorrelation analysis in [77]. (2) Gaussian process (GP) based non-parametric time-series prediction in [78]. (3) Maximum likelihood classification in [100]. (4) Joint fully convolutional auto-encoders (FC AE) model in [66]. The results suggest that our method provides significant higher levels of both accuracy and recall in comparison to the methods based on SARIMA, GP, and AE. Unfortunately, [100] did not report values in regard to recall.

Focusing the discussion specifically on the study presented by Sublime and Kalinicheva in [66], they proposed a deep-learning method for change detection (FC AE). Their method is useful for post-disaster damage mapping. Although they also deal with real-world problems, their study differs from ours because their method detects outliers rather than apply the taxonomy [13] (i.e., detecting anomalies) to solve real-world problems. According to Sublime and Kalinicheva, their outlier detection method achieved the best results when it was compared against other machine-learning methods. The levels of accuracy and recall achieved by our strategy are significantly higher than the ones achieved by the deep-learning method proposed in [66]. Additionally, the level of precision achieved by the deep-learning method proposed in [66] is lower than the one achieved by our strategy. Therefore, the results suggest that our method provide results that are more efficient and relevant, and a higher quantity of truly relevant results compared to machine-learning methods applied for post-disaster damage mapping, for example.

As previously mentioned, to the knowledge of the authors this is the first endeavor of a series of studies dedicated to investigating the detection of anomalies in accordance with the taxonomy published in [13] and their potential applications to remote sensing. This study proposed a strategy taking into consideration the multifaceted nature of anomalies [13], whereas other detection methods have been dealing with different types of anomalies as if all of them were outliers. Since the present strategy is based on this taxonomy [13] which presents ten different domain anomalies, this study overcomes the number of categories studied in previous investigations, e.g., [25] which presented four different anomalies. The smallest differences regarding the recall in comparison to our study were reported by Qi Liu et al. in [25] and Chen et al. in [82]. However, Chen et al. did not report their achieved accuracy. Overall, our study provides relevant improvements over other studies and also important support and guidelines for future studies.

Regarding limitations, it should be noted that the accuracy of the proposed strategy varies depending on the environment of application. It is inevitable when many different geographical features are evaluated by some anomaly detection strategy. However, the differences in levels of accuracy of the proposed strategy are of little significance. Further work is planned to examine the effects of applying filters on the images before applying them to the proposed strategy to determine if the different levels of accuracy can be minimized among different environments. The results indicate that the proposed strategy accurately detects anomalies related to brown mud and crude oil in the majority of evaluated images from remote sensing.

Regarding negative results, it should also be noted that the detection of incongruence, when the water pollution is caused by either red mud or sewage (Table 6 and Table 9), was not obtained by the proposed strategy. This limitation is probably caused by the use of seven bands as a pattern for all experiments of this study. However, it is of little significance, since the ability of the proposed strategy detecting incongruence was demonstrated when the water pollution is caused by either brown mud or crude oil. All in all, the achieved results indicate that the proposed strategy detects the presence of anomalies of the type unexpected structure and structural components, in accordance with [13].

8. Conclusions

This paper presented an incongruence-based anomaly detection strategy for analyzing images from remote sensing, with the aim of describing its practical application in detecting water pollution. Therefore, the practical application of incongruence was introduced as a strategy for detecting anomalies in real-world applications, i.e., non-synthetic data present in images from remote sensing.

The proposed strategy detected anomalies of the type unexpected structure and structural components, distinguishing them from outliers. All processing was performed on just a single Landsat 8 image, whereas other methods use time-series for each environment to be able to detect outliers. This means that systems dedicated to observing or monitoring occurrences of anomalies of this type will save time and use few computational resources because they will process a single image. Therefore, systems based on the application of our strategy will help by saving more money than systems based on other proposed detection methods.

This study has established a relationship between the occurrence of anomalies of this type and the presence of brown mud or crude oil as either potential contaminators or contaminators of large rivers. Our strategy identified incongruences to detect oil spills in some areas where there is the presence of water, or brown mud where there are high levels of turbidity in the water. These findings are relevant since an incongruence-based anomaly detection strategy can significantly increase the ability of surveillance systems for detecting environmental disasters and for performing mappings.

Although the achieved results are encouraging, future studies should investigate: (1) The effects of using Gaussian process-based filters on the images before applying them to the proposed strategy in order to know if it would be useful to minimize differences between the levels of the accuracy among different environments. (2) The effects of applying feature extractors on the images as a preprocessing stage before applying the proposed strategy in order to determine if the classification process can be improved for all environments. (3) Whether our strategy will be able to detect other pollutants not using pan-sharpening.

Future work will investigate: (1) Practical applications of other types of anomalies individually, such as unknown object, measurement model drift, unknown structure, unexpected structural component, and component model drift. (2) A single system to deal with the practical application of all types of anomalies together.

It is expected that this study will open up an entirely new range of applications for computational tools on images from remote sensing in order to model and analyze environments. Moreover, this anomaly detection strategy can be applied to a wide range of studies in geosciences and other scientific areas.