**Smart and Climate-Smart Agricultural Trends as Core Aspects of Smart Village Functions**

## **Adegbite Adesipo 1, Oluwaseun Fadeyi 2,3, Kamil Kuca 3, Ondrej Krejcar 3,4,\*, Petra Maresova 5, Ali Selamat 3,4 and Mayowa Adenola <sup>6</sup>**


Received: 1 September 2020; Accepted: 17 October 2020; Published: 22 October 2020

**Abstract:** Attention has shifted to the development of villages in Europe and other parts of the world with the goal of combating rural–urban migration, and moving toward self-sufficiency in rural areas. This situation has birthed the smart village idea. Smart village initiatives such as those of the European Union is motivating global efforts aimed at improving the live and livelihood of rural dwellers. These initiatives are focused on improving agricultural productivity, among other things, since most of the food we eat are grown in rural areas around the world. Nevertheless, a major challenge faced by proponents of the smart village concept is how to provide a framework for the development of the term, so that this development is tailored towards sustainability. The current work examines the level of progress of climate smart agriculture, and tries to borrow from its ideals, to develop a framework for smart village development. Given the advances in technology, agricultural development that encompasses reduction of farming losses, optimization of agricultural processes for increased yield, as well as prevention, monitoring, and early detection of plant and animal diseases, has now embraced varieties of smart sensor technologies. The implication is that the studies and results generated around the concept of climate smart agriculture can be adopted in planning of villages, and transforming them into smart villages. Hence, we argue that for effective development of the smart village framework, smart agricultural techniques must be prioritized, viz-a-viz other developmental practicalities.

**Keywords:** smart village; smart agriculture; climate-smart agriculture; technology; sustainability

## **1. Introduction**

The need to develop rural communities in terms of productivity and convenience, so as to curb urban migration has received much attention in the last decade. First, the Institute of Electrical and Electronics Engineers (IEEE), as part of its mission, commenced the installation of solar-powered bulbs in many rural communities worldwide [1]. This was followed in 2016 by the Cork Declaration, agreed amongst 340 representatives of European states towards ensuring that rural communities enjoy better lives. These efforts culminated into the coining of the word "smart village", defined as a

community that tries to develop current strength and resources, while making futuristic developmental plans on the basis of technology [2,3]. While there are several thematic areas of priority within the smart village development framework, agriculture is seen as the most important of them all [3]. Furthermore, the need to bridge the digitization gap between cities and villages, is also an important aspect, so that lives and livelihood can be improved. Since a smart village is one that seemingly accepts new technologies, precision agriculture uses ultra-modern techniques for animal and crop production, which saves time and reduces wastage, and meets the requirements of smart villages. This is crucial for the sustainability of smart villages [4]. This is because improved food production and efficient animal management systems must be at par with village development, and must be continually transformed to influence the different aspects of smart villages, in terms of policy and practice [5].

To effectively play its role in smart villages, precision agriculture covers smart and climate smart agriculture (CSA) techniques, and other aspects that are capable of ensuring higher agricultural production output in an environment-friendly manner, provides optimum income for the farmer, and is able to feed a growing population. Many studies showed that these processes can be realized through the adoption of ultra-modern agricultural techniques such as bio and nano technologies [6], IoT and blockchain-based methods [7], and drone technologies [8], among other climate smart ideas. On the basis of this argument, efforts that tend to reduce farming losses, increase yield, as well as monitor, detect, and potentially prevent plant and animal diseases are now being automated, finding growing applications, and offering optimal solutions. Based on the forgone explanations, the current study attempts to establish smart and CSA trends in smart village research, in order to see how much they are useful for smart village development.

The rest of this study is arranged as follows. Section 1.1 draws a foundation for this study, by focusing on the research question. Section 2 briefly builds a background for smart village research by listing existing projects, and describes a few state-of-the-art smart agricultural solutions. In Section 3, attention is drawn to climate-smart agriculture, with specific reference to what makes up the concept, a few challenges in its framework, as well as the latest progress in its development. Section 4 describes the challenges created by the interplay of adopting CSA in smart villages, and also tries to answer the research question. The section also conceptualizes climate-smartness, as it influences sustainable development of smart villages. Finally, Section 5 describes future research directions in smart-village and smart agricultural research, and draws relevant policy recommendations and conclusions

#### *1.1. Research Question*

Based on the vast importance of agriculture in smart village development, this study adapts its research question from the editorial note presented by the editors of MDPI's special issue within the Sustainability journal published in August 2018. Within the report, Visvizi and Lytras [9] gave a revealing background of future directions for smart village research. The editors pointed to a few research questions that future smart village research should strive to answer. One of these is: "How will smart and CSA research give account of, and conceptualize transformation and change in the smart village context?"(p. 8) [9]. This question is what the current study modifies and seeks to answer.

## **2. Related Literature**

## *2.1. Current Smart-Village Projects Around the World*

Before delving completely into smart agricultural systems in smart villages, it is important to consider existing smart village initiatives in order to have an updated knowledge of smart village trends, and why smart agriculture might need to be prioritized. Zavratnik et al. [10] described the IEEE smart village project, and the EU smart village initiative, which are further elaborated in subsequent paragraphs.

The IEEE smart village program is one of the most popular today. It has a goal of advancing education in off-grid societies, and fostering sustainability in the entire value chain of the smart village energy sector. Initially taking off as an initiative that seeks to provide community solutions in 2009, the current name was coined 5 years later. The IEEE smart village plan is a global initiative, touching lives in Asia, some parts of North America, and mostly in Africa [1], through the promotion of smart energy production in rural areas, and is mostly financed through fundraising. Major efforts that were developed from the initiative include the so-called SunBlazer II—a movable power base solar station [11]; "Learning beyond the Light Bulb" [12]—a program aimed at training locals on the development and design of off-grid solar electricity panels and fostering its sustainability and scalability. As reported by Zavratnik [10], the program also comes with a remote study event that runs for about nine months, and allows practice exchange amongst involved communities, for knowledge sharing and skill enhancement.

Within the framework of the Consultative Group on International Agricultural Research (CGIAR), several smart-village projects took off around the world. Many of which were funded through international research organizations with clear impacts in areas that were worse hit by climate change [10]. These projects mostly focus on training smallholder farmers on agricultural resilience, through the adoption of practices that support food security [13], so that persons within these affected communities are able to maintain a livelihood through agricultural methods that help decrease GHG emissions. For instance, farmers in the Lower Nyando valleys of Kenya are benefitting from improved agroforestry systems that adopts knowledge of Information and Communication Technology (ICT) [13]. Thanks to the CGAIR initiative, they are able to cultivate cash crops in-between rows of multi-purpose trees, thereby improving soil stability and enrichment. Given the increased demand for trees, several nurseries were developed, adding farmers' incomes, and with women as the highest beneficiaries. The state of Bihar in India also benefitted from the CGIAR's smart village initiative. It previously had soils that were greatly affected by water-logging, but new drainage construction changed the channel of rapidly flowing flood waters, out of the farming areas [14]. This improved system also ensured that underground aquifers were steadily recharged. Improved technological ideas also saw better rainwater harvesting in areas benefitting from the Climate Change Agriculture and Food Security (CCAFS) program. Overall, weather and planting can now be monitored from smartphone applications by the farmers, in order to avoid unwanted losses [15].

The European Union's smart village initiative is by far the most organized and detailed system. Having undergone several fine-tuning, the initiative has improved tremendously since the Cork Declaration 2016 [10]. Notable amongst the goals of the EU smart village drive is agricultural boost, mainly because the rural areas are where European foods are mostly produced [16]. There is also the goal of reduced youth exodus to urban centers [17]. In the first assembly of the newly adopted "Intergroup SMART Villages for Rural Communities", György Mudri, a former Members of European Parliament stressed that smart villages are not only for the development of new infrastructures, but also for building capacity of locals [18]. In response to this statement, The Austrian Chamber of Agriculture commenced online training for about 10,000 farmers, who now have remote access to latest agricultural researches and can subsequently implement such ideas on their farms [19]. There is also the so-called COWOCAT rural initiative, which currently trains youth to commence working in villages [19].

Description of smart village drives of the above initiatives show that agriculture is one of the most prominent aspects of the smart village plans. As a result, this study delves into smart agricultural practices that can the build capacity of smart villages, if adopted.

#### *2.2. Ultra-Modern Smart Agricultural Solutions*

While there are varieties of smart technologies adopted in agriculture nowadays, this section focuses only on bio-sensors, agricultural drones, IoT and Blockchain-based sensors, and a number of combined technologies that adopted animal husbandry, as well as in crop, soil, and pest management.

#### *2.3. Nanostructured Biological Sensors*

Biological sensing devices are some of the new technological interventions reshaping agricultural systems today, which might be adopted in smart villages. As reported by Antonacci [6], bio-sensors with extremely small structures were found to possess the ability to help in crop maturity evaluation, management of amount of pesticides and fertilizers, as well as detection of humidity levels in soils for effective irrigation. To carry out these functions, bio-sensors rely on the characteristics of nano-materials, such as immobilizing bio-receptors on transducers, integrating and miniaturizing some biological components of plants, transducer systems, and micro-fluids, into very complex plants make-up [20,21]. Although the use of harmful pesticides is gradually being phased out in agricultural systems, less harmful pesticides are still very much in use [22]. In areas known for prior pesticides usage, modern agricultural techniques often aim at detecting pesticide presence, as well as their levels within the soil, before cultivation. To do this, cutting-edge bio-sensors with very high sensitivity (because of their surface-volume ratio), extremely rapid response time, and quick electron-transfer kinetic are utilized. The sensors possess stable strength to map pesticide quantities within soils, and longer lifespan, when compared to the earliest bio-based sensors [23]. Newly improved bio-sensors with extremely small structures are also able to surpass soil pre-treatment, due to the presence of pesticides, herbicides, and fungicides, without losing their potency [6].

Yu et al. [24] developed tyrosinase/TiO2 biosensor to determine the presence of atrazine pesticides. This was done by fabricating a structure through the allowance of vertical growth of TiO2 nanotubes. This meant that well-arranged nanotubes would provide large surface areas for immobilizing the tyrosinase enzyme. The structure gave room for excellent loading of enzymes, as well as transfer of electrons, which yielded improved system robustness and sensitivity. The system was tested in well-grinded, air-dried paddy soils, gathered at varying depths. The soil also passed through a sieving process using a 1.0 mm filter, and a 35 ◦C re-drying process that lasted for 48 hrs. It was subsequently mixed with acetone, prior to undergoing shaking at a temperature of 25 ◦C for 60 min. Results given by [24] showed that after carrying out analysis of supernatants, atrazine was observed to be present in 0.2 ppt to 2 part-per-billion. Standard deviation was subsequently found to be below 0.05 ppt when compared to high performance liquid chromatography (HPLC).

Dong et al. [25] introduced a novel nano-structured bio-sensor technique for detecting very low pesticide traces in soil. The technique works by electrochemically reducing Ellman's reagent via the inhibition of acetylcholinesterase. This bio-sensor adopts amperometric, designed to immobilize acetylcholinesterase on multiple walls of carbon-type nanotubes-chitosan nanocomposites modified glassy carbon electrode. High sensitivity of the system is offered by the very good conductivity and biological compatibility of multiple walls of carbon-type nanotubes-chitosan [25]. This can be additionally improved by electrochemically reducing 5,5-dithiobis (2-nitrobenzoic) acid. In testing the system, methyl parathion pesticide was observed to exhibit an inhibitive effect on acetylcholinesterase. An electrochemical change in the reduction response of 5,5-dithiobis (2-nitrobenzoic) acid was also observed. Overall, the system was found to possess a pesticide detection precision of 7.5 <sup>×</sup> 10−<sup>13</sup> M when tested on spiked soil.

In another nano structured bio-sensor study, Shi et al. [26] observed the presence of soil acetamiprid using SELEX; a new 20 mer bio-sensing unit that is able to bind acetamiprid, using aptamer made of nanoparticles of gold. The unit works to detect the pesticide optically at values ranging from 75 nM to 7.5 μM. It bears the combined characteristics of a nanomaterial, and those of artificial molecules. Tested soils were collected around Tongji University, China, with initial air-drying carried out before the sample was grinded, to allow 1.0 mm sieving. A second drying was also done using an oven at 35 ◦C for 2 days, prior to acetone mixing and shaking at 25 ◦C for 60 min. Dichloromethane was subsequently added to the mixture, and then removed ultrasonically before the sample was filtered.

Beyond sensing pesticides within soils, bio-sensors with very tiny structures were also employed in monitoring diseases of crop plants. A very important aspect of any smart village is the effective management of farm economy, achievable through the protection of crops against diseases. Quantum dots offer classical examples of materials that are useful for monitoring plant diseases, as they possess broad excitation spectra. Safarpour et al. [27] identified the vector responsible for sugar beet's yellow vein and Rhizomania disease like Polymyxa betae. This was detected using quantum dot techniques that subjected the plant root sap samples to several pre-treatment in order to extract the virus. The quantum dots unit utilizes Förster Resonance Energy Transfer (FRET) modeling in its detection operation [27]. By using a similar technology, Bakhori et al. [28] detected synthetic oligonucleotide of Ganoderma boninense. However, this work employed adjusted quantum dots with carboxylic groups that are then conjugated using a DNA probe. This gave rise to an improved sensitivity of the system, yielding 3.55 <sup>×</sup> <sup>10</sup>−<sup>9</sup> M as the detection limit [28].

By adopting bio-sensors in the detection of soil nutrients and fertilizers, Ali et al. [29] revealed that soil nitrates can be detected using a system that relies on microfluidic impedimetric sensing. The unit works by adopting nano-sheets of graphene oxide and the nanofibers of the so-called poly (3,4-ethylenedioxythiophene). The researchers showed that poly (3,4-ethylenedioxythiophene) composite can bear the enzyme; nitrate reductase, and also measure the amount of nitrate ions in soil samples on which sweet corn was cultivated. This is done at 0.44-442 mg/L concentration, so the detection limit was 0.135 mg/L [29]. Carrying out the procedure, however, involves sample drying at 105 ◦C, and subsequent nitrate extraction through the addition of 2 M KCl solution. The mixture was shaken for 60 min, and filtered using Whatman filter paper. Finally, sample extraction was kept in a syringe for infusion into the experimental device.

Several other research examples exist for bio-sensor utilization in agricultural work. Nevertheless, a summary of some state-of-the art techniques, some of which are already described elsewhere, is presented in Table 1.



#### *Sensors* **2020**, *20*, 5977

### *2.4. Drone Technologies (Unmanned Aerial Vehicles)*

Unmanned Aerial Vehicles (UAV), also known as drones, have become popular in agricultural production work. In a review study by Mogili [8], the researchers reported that drones can be used in pesticide and fertilizer application, so that humans do not come in contact with the some of these pesticides, which are harmful, and are gradually being phased out. Drones can also function as water sprinkling systems [8,36].

Primicerio et al. [37] adopted VIPtero, a UAV for managing a vineyard in an experimental set-up in Italy. The system, which is made up of an aerial platform with six rotors and a camera, can fly in a self-governed manner to a particular point in the air, in order to take measurement of the vegetation canopy reflectance. Prior to flight take-off, accuracy of the camera is evaluated in relation to ground-based measurements with high resolution, which were gathered using field spectrometer. Subsequently, VIPtero gets air-bound in the vineyard, and gathers as many as 63 multi-spectral images in a 600 seconds time period. The recorded images are analyzed and classified, prior to the production of vigor maps on normalized difference vegetation index. Results showed the heterogeneity conditions of the crops, implying that they were in line with those gathered using the ground-based spectrometer [37]. This smart system appears to be promising as an effective and detailed data gathering system in agriculture, and can be adopted over larger areas in smart villages.

In another UAV based research, Burgos et al. [38] used a 4 cm Sensefly Swinglet UAV to differentiate green cover from grape canopy. A digital surface model (DSM) with 3 dimensions was adopted to create an exact digital terrain models (DTM), acquired via the use of processing libraries of python, and subsequently subtracted from DSM, so as to arrive at a differential digital model (DDM) for the measured terrain (a vineyard). Vine pixels within the DDM were obtained by selection of pixels >50 cm elevation from the ground. The results indicated that there is a possibility of separating vine row pixels from green cover pixels, as a differential digital model pointed to values ranging from –0.1 m to +1.5 m. Furthermore, manual polygon delineation, which depended on an RGB image of the vine rows and green cover, revealed huge differences averaging 1.23 m and 0.08 m for vine and ground, respectively. Elevation of the vine rows was good and tallied with its topping height of 1.35 m from the field [38]. The authors noted that vine pixels extraction would aid future analyses, such as pixels' supervised classification.

Berni et al. [39] also demonstrated the possibility of generating remotely sensed data over an agricultural field, using a UAV that had a relatively cheap narrowband and thermal multispectral imaging sensors of 20 cm and 40 cm resolutions, respectively. The system gave rise to surface reflectance and temperature data, after adapting MODTRAN-based atmospheric correction. Biophysical parameter estimation was carried out using a number of vegetation indices, leading to the production and validation of chlorophyll content, detection of water stress from PRI index, as well as the temperature of the canopy. These results showed that the system yielded the same results as the conventional, expensive, and risky manned airborne sensors.

## *2.5. IoT-based Sensors with Complimentary Blockchain Technology*

Many villages face severe agricultural challenges and that require upgrading to smart agriculture, which offers a wide range of state-of-the-art solutions. For example, in villages where access to water is a challenge, Khoa [7] maintained that IoT-based sensors can be useful in water-management irrigation systems on large rural farms. In their research, the authors developed a novel system that is able to monitor soil water level and schedule sprinkling/spraying times in well-calculated amounts. This relatively cheap technique, functions by receiving real-time data from sensors fixed within strategically arranged tunnels, in and around the farm. Based on the information supplied by sensors, which can be received through a mobile phone application, the user might decide to water the farm. Subsequently, when soil water level increases to an optimal level, the system notifies the user, who can remotely or manually switch-off the water-pumps. A unique feature of this system is its usability in up to two farms [7]. In a similar study, Nagpure et al. [40] described another IoT-based

system that works by using a similar routine as [7]. However, two differences include; scaring animals away using current pulses, and wireless sensor monitoring of the ecological conditions (e.g., altitude and humidity) to ascertain the amount of irrigation water needed each time [40], which the latter unit possesses.

Mat et al. [41] presented an IoT-based mushroom cultivation, which produced a better yield when compared to the conventional system. This tool is based on an automated sensors for fertilizer application and water sprinkling on the farm, which can be controlled from the farmer's mobile phone or manually, from a centralized point within the farm. The system ensures that the timing for wetting the crop is strictly adhered to, so that the farming operation can progress even without the farmer. Overall, it was observed that average mushroom size in thickness and weight exceeded conventional cultivation by 0.3 cm and 5 gram, respectively.

As reported by Prathibha [42], it is important to curb the effect of environmental conditions on crop yield output. To do this, an efficient measurement method of the elements of weather might be required. Prathibha's research, therefore, proposed a CC3200 combined sensor unit, which comprises a processor for network, a micro-controller, a Wi-Fi unit, a camera, as well as temperature and humidity sensors. This weather utility device comes as a portable unit with low power consumption for longer battery-life. The system monitors temperature (using a thermopile sensor that uses infrared technology) and humidity across the agricultural field, which are subsequently processed as camera images and sent via Wi-Fi to the farmer's mobile phone as multi-media messages. Information of this nature helps the farmer to know how good the soil water is to support the grown crop. A similar study [43] designed another unit that can also send immediate signals to a farmer, after recording real-time data on weather, in and around the farm. This unit was made up of a breadboard, a combination of sensors that can monitor UV Index/IR/Visible light (SI1145 Digital Sensor), soil moisture content, humidity, temperature (DHT11), an ESP32s Node MCU, all of which are connected to a monitor, which is in turn linked remotely to the famer's mobile phone, using an LED visual alert and Blynk mobile phone application. Two very special characteristics of this unit are; its ability to save power in sleep mode, for a battery life that averages 10 days, and the speed of sending signal (180 s) [43].

By using a pre-coded algorithm, also known as the "*Cuckoo Search Algorithm*" [44], a framework for automated watering of a piece of farmland was designed. Based on pre-analysis of different kinds of soils, the researchers found that a soil moisture value of 700 meant dry, and would require immediate watering. The IoT sensor was, therefore, designed on the basis of this information, comprising the so-called ThingSpeak, which also gives direction on the most suitable soil type for a specific crop. A temperature-based sensor was initially used on the soil, and the result was sent to a converter called Arduino. Depending on the measured value, the Arduino was connected to an automated watering system, which could be controlled from the mobile device handled by the farmer. When soil wetting gets to optimal levels, the soil sensor sends a feedback signal to the farmer who then stops the watering [44].

While IoT-based sensing techniques are available for improving precision agriculture, optimizing these ideas with blockchain technologies might offer even more robust results. Patil et al. [45] noted that IoT-based sensing technologies might sometimes be flawed on the grounds of; extremely large scale, a lack of homogeneity of different IoT-based sensing operations, as well as standardization. This meant that the data gathered using IoT-based sensing technologies comes with privacy concerns for the farmer [45]. Hence, the researchers developed a blockchain greenhouse farming tool to cater to security and privacy. The model is made up of a smart greenhouse (a covered piece of farmland protected from environmental conditions) that comprises a series of sensors and actuators, smart hub (a local blockchain which manages the connectivity of all sensors and equipment in the smart greenhouse); an overlay network connection that manages the nodes; a cloud storage platform and an end-user platform. A system of this nature addresses security challenges across all fronts within the farm [45]. In another Blockchain-based smart agricultural study [46], a traceability platform for food safety was designed. In collaboration with the Internet of Things, the system involved "Enterprise Resource

Planning" where farmers, processing plants, and organizations involved in the logistics of agricultural and food products, and the consumers, can assess on their mobile phones, a blockchain node that gives detailed description of how the products were cultivated, harvested, stored, processed, and sold [46]. The essence of this technique was to build virtual trust in food processing, using the so-called *Trusted Trade Blockchain Network Cloud Platform (TTBNCP)*

#### *2.6. Smart Animal Production, Management, and Monitoring*

The use of machine vision in body condition scoring of dairy received extensive research attention in the last few decades. Fox et al. [47] listed animal nutrition, insemination, and health as core reasons for body condition monitoring. When this process is carried out by the farmer or veterinary doctor (human monitoring), there is a possibility for biased data gathering, due to the individual's mental state, level of experience, and residual knowledge [48]. Furthermore, the process might also be time-consuming. Improvement in body condition scoring in the 1970s employed ultrasounds [49], which was flawed on the ground that mastering the collection of reliable ultrasonic body condition scoring required more time in comparison to data gathering by humans [50]. Additionally, the cost of purchasing the ultrasound device, and hiring an expert, made the process too expensive. This led to camera-recording of animals applied to body condition monitoring, based on the belief that this would yield better results [48]. Fourier descriptor cameras [51], thermal [52], and RGB cameras [53] were adopted. Nevertheless, images could not be processed automatically until seven years ago [54]. As a step in the right direction for body condition scoring, Spoliansky et al. [55] used a 3D camera that was well-equipped to carry out automatic image processing, leading to the development of effective and unbiased collection of body condition scoring, in 2017 [55]. The system provides real-time data, useful for commercial milking purposes, and genetic evaluation (based on lactation). In addition, automated image gathering of body condition scoring might provide ease of monitoring when there are more than one animal. This is pivotal to early warning signs for morphological changes in animal body size [48].

Yanmaz et al. [56] suggested adopting thermography in the early detection of lameness in horses. Two types of thermography are—contact and contactless types. Contactless thermography has higher precisions via infrared radiation. Internal temperature of the affected animal part can be viewed on a medical thermogram, so that treatment can be planned early. Similarly, temperature around a sick cattle's gluteal region differs significantly from those of other parts when studied using thermal infrared scanning [57]. As reported by Steensels et al. [58,59], temperature management in poultry, as well as early mastitis identification are some other areas where thermal scanning was reported to be useful.

Accelerometers exist nowadays for remote measurement of animal gait [60]. This can determine when the animal is lying or walking [48]. The method is also useful in the determination of lameness, animal sickness, or how much the animal feeds [61], based on the distance covered by the animal against time. The device can be mounted on the leg, ear or other parts of the animal, so that it continuously sends signals to the farmers' mobile phone. When the animal is lying down, the instrument automatically changes its processing speed. The device was also utilized in monitoring the health of fishes, by attaching it to their fins [62]. A demerit of the system, however, is the fact that continuous processing of mobility data tends to rapidly reduce battery life [63].

Wearable belts that are able to tap animal sweat and measure the amount of sodium it contains are some of the latest smart technologies in animal husbandry. In the work of Glennon et al. [64], the authors developed a smart technique for quick detection of the sodium content of sweat. The unit appears in two ways, one that resembles a vertically placed watch and can be worn round an animal leg, and the other that looks like a horizontal pod. Both come with a Velcro strap that can be used to attach it to the animal skin. The systems, through capillary action, receive sweat via its orifice, and send it through an electrode that is sensitive to sodium. The electrode in turn sends it into a storage section containing an adsorbent substance. Sweat flow rate can be improved by varying the width of the sweat

flow channel in-between electrode and the storage section of the system. Sweat flow rate generally decreases with decreased width of the flow channel. This also determines the length of time the system will be used before the electrode and the adsorbent material are changed. Stored sweat is available for measurement as total harvested sweat volume as well as its sodium concentration per time. Electrode signals are moved to an electronic board that possesses high input impedance for voltage capturing. Results of this analysis are sent to a remote base station, which is either a laptop, or mobile phone, using bluetooth technology for onward visualization and possible storage.

In livestock management [65] a pregnancy detection method was developed, which makes use of Xbee transmitters linked to LM35 temperature sensors. Two animals were experimented, with temperatures recorded at five days and twelve days from insemination. The sensing system was attached to the tail of the cows. Temperature records were found to be high in pregnant cow, especially in the evenings. The sensing unit works effectively within a distance of 40 m, and serves as a low-cost technique, when compared to some invasive pregnancy detection methods in livestock.

In an ongoing study aimed at finding a sensor that is able to detect the level of progesterone in milk fed to cattle [66], interdigital sensors are used, which are able to yield single side access to the substance being tested. Within the study, progesterone hormone of about 20 mg was allowed to dissolve in 0.5 mL of approximately 100 per cent ethanol. The solution was subsequently poured into about 1000 mL of pure Milli-Q water, so that a stock solution was achieved. Successive mixture dilution led to 0.02 ng/mL progesterone concentration. The sensor offered different results at varying progesterone concentration with sensitivity in the pitch range of 50 μm [66].

Infectious coughs in piggery need rapid detection and treatment. To detect this kind of cough, Ferrari et al. [67] fixed 1 m multi-directional microphones of 50 to 16,000 Hz around a farm. The microphone was connected to a laptop, and animal cough sound patterns were recorded, digitized, and analyzed using Matlab 7. Having earlier injected healthy animals with citric acid, acoustic parameters such as time difference between coughs, peak frequency, and root mean square were used to differentiate coughs from a healthy pig from those of infected ones. While healthy pigs relaxed for about 52 s after each cough attack, the infected pigs coughed after every 37 s. Peak frequency for infected and non-infected pigs was observed to be 1600 Hz and 600 Hz, respectively.

### **3. Moving towards Climate-Smart Agriculture**

Having established in Section 2.1 that agriculture is one of the most important factors to be considered in smart village development, it is crucial to stress that climate change is a major stressor for agricultural development of rural communities [68]. The implication is that developing an agriculturally-smart village entails accepting the concept of climate-smart agriculture. Agricultural risk posed by climate is a threat to food security. As a result, there is an urgent need to effectively manage agricultural production, while fighting climate change through adaptation, resilience and mitigation [69]. This is what climate-smart agriculture offers.

There is currently no unified definition for climate-smart agriculture (CSA). In fact, almost every new study within the framework of smart-agriculture, views CSA in a slightly unique way. Nevertheless, to build a strong foundation for climate-smart agricultural framework in smart village development, the current study adopts existing knowledge and definitions, to coin a new and more robust definition for the term. Table 2 presents some definitions put forward by climate change and agricultural scholars and research organizations. Keywords derived from the definitions show that each has one or more shortcomings. As a result, it might be difficult to build the concept of smart village on a definition that lacks one or more fundamental aspects.

Given the definitions in Table 2, considerable aspects of climate-smart agriculture include; capacity building, sustainability, emission reduction, vulnerability reduction, profit, food security, transformation, new knowledge, new technology, and productivity. By linking the above keywords together, we define climate-smart agriculture as a "transformative and sustainable kind of agriculture that tries to increase efficiency (productivity) in food security and production systems, using a

combination of the pillars of climate change (adaptation, resilience, and mitigation) as well as smart and new technological knowledge, that do not only build capacity of farmers' in terms of farming techniques, but also increase profit, reduces vulnerability of the systems as well as their results (farm products/animals), through the reduction of GHG emissions."


## **Table 2.** Definitions of CSA.

While it can be argued that the list of keywords suggested within the current study is not exhaustive, many other definitions tend to be built around at least one of these keywords. Figure 1 is a diagrammatical representation of the main aspects of climate-smart agriculture for which it stands as a significant part of a smart village. The implication of the above expository listing of the fundamental parts of climate-smart agriculture means that for a smart village to be so called, it must strive to maintain within its agricultural systems all different aspects of CSA. Furthermore, other aspects of the smart culture within the smart village setting; smart energy management, smart living and smart healthcare, etc., must tap from these fundamental attributes of CSA, in order to provide robust services in their smart village functions.

**Figure 1.** Key aspects of climate-smart agriculture (CSA).

In demonstrating whether CSA could increase rice yield in China, Xiong et al. [79] used crop simulation models; version 0810 of the Environmental Policy Integrated Climate (EPIC) model [80], and version 4.0 of the so-called DSSAT, an acronym for Decision Support System for Agro-technology Transfer [81], respectively. It was observed that these software simulations that gave ideas on cultivar improvement and optimization of management practices for rice due to climate change, led to increased rice production. The EPIC models specifically yielded over 2000 kgha−<sup>1</sup> during the 30-year period under review [79].

Rural African farmers tend to suffer a lot from adverse weather conditions. This further creates a need for cheap and reliable weather forecast system. To attend to such needs in Nigeria [82], a cheap automatic weather station that functions on solar energy was designed. By linking meteorological sensors to microcontrollers, the farmer could gain access to processed information related to weather, through a television screen. A thermometer collects temperature information, while the anemometer and LDR measures wind speed and sunlight, respectively. Embedded temperature sensors within the microcontroller receives analog information gathered by the thermometer and converts it to digital signals [82]. In some cases, unprocessed data can also be sent to farmer's mobile phones. The cheap rate of the unit shows that it can serve as a very good system for crop management and food security, in the least developed nations.

In a research carried out by Tenzin et al. [83], to ensure effective weather monitoring around a farm, the authors designed a very cheap cloud-based weather measurement unit, using an integration of different unique weather sensors. The system, which is made up of a base and a weather station, as well as a display unit, is capable of effectively gathering humidity, temperature, wind direction, wind speed, and many other weather data types. By experimenting its usage and statistically analyzing gathered data, it was observed that the unit provided similar results as the Davis Vantage Pro2 weather monitor, which was pre-installed on the same farm, thus, offering a cheaper option [83].

In a bid to design an integrated farm that efficiently manages water and reduces climate-demanding inputs, Doyle et al. [84] designed an aquaponics unit for vegetables and fish. The design consists of a 12V DC pump that delivers water from the fish tank to the flood tank, which then supplies the area where the crops are planted at a constant rate. As soon as water is removed from the fish pond, it is carried by gravity through the grow bed area, where it is stored until it is needed for watering the vegetable bed. The pump is powered using a solar panel of 150-Watt with a 120 Ah battery.

Having described some smart agricultural and climate-smart agricultural studies, it is important to note that while smart agriculture is mostly developed, research on CSA is relatively new and still at the level of policy and framework description [85]. In a systematic review study by Chandra et al. [85], the authors observed that research on CSA is mainly divided into three parts; global policy and plans around the world concerning further development of the concept, scientific research directions, and integration of pillars of the concept (which includes; adaptation, resilience, mitigation, and food security). With respect to CSA policy framework developed by the World Bank, Taylor [86] faulted the fundamental make-up of the concept on the following grounds.


Given these fundamental shortcomings of CSA [86] 'climate-wise food system' is suggested as a more direct term that should be used to refer to sustainable food production systems, rather than CSA. Another criticism on the policy and framework of CSA comes with the injustice meted to smallholder farmers, as a result of the implementation of the concept [87]. By administering interview to some CSA experts, analysis based on a number of ethical positions showed that implementation of climate-smart agricultural approaches is not fair, especially with respect to allocation of income benefits and challenges of cost associated with emission reduction [87], among smallholders farmers and small agricultural processing industries. Budiman [87] further argued that based on how climate justice works, sharing of income benefits should depend on the financial capability of farmers.

In a comparative study of Philippines and Timor-Leste, five important features of climate-smart agricultural practices were observed by Chandra and McNamara [88]; strategies at country-specific institutional levels; delegated financial procedures; the state of the market; technology; and knowledge. In the two countries, CSA was used to resolve climate vulnerability challenges more than it was associated with emission reduction goals [88]. Overall, the researchers observed that advancing the course of CSA in these countries might involve multi-stakeholder approaches that cuts across different levels of participation, both within and outside the farm, rather than mere technical CSA developmental inputs [88]. From the above arguments for and against CSA, it is clear that while there are still fundamental challenges revolving round the CSA concept, the terms might likely continue to be utilized for agricultural problem solving, until it attains uniformity and intersection of ideologies, amongst researchers and policy makers.

## *What does Smart- and CSA O*ff*er Smart Villages?*

Having described in previous sections how the concept of CSA has evolved amidst the challenges faced within its developmental framework, an examination of the utility of climateand technology-driven agriculture to smart villages is important. According to Azevedo [5], there is a big chance that CSA will empower and strengthen the conceptualization and execution of smart village in different ways. Safdar and Heap [89] noted that development of small grids to power certain climate-smart technologies has so far spurred a re-imagination of the possibility of home solar powering in many Indian villages. Items such as solar lanterns, and street solar lighting systems have become very popular. Nevertheless, a new concern is the way to enhance local productions and repairs of these materials, in order to cater for higher tariffs of importing them to interior villages, and shipping them

back for repairs, when the items develop technical faults. The report also stressed how CSA has so far upheld gender equality, for instance, the CCAFS project in Kenya's Nyando valley has mostly favored women whose incomes have improved due to new technology for growing their vegetables [13].

In documenting how CSA could provide smart village farmers with possible economic benefits, Khatri-Chhetri et al. [90] carried out a research using farmers of India's Indo-Gangetic Plains. Major CSA practices by the farmers include diversifying crops, land levelling using laser, nutrient management in a site-specific mannerism, management of residue, and zero tillage, among others. The researchers started by calculating how much the farmers spent to adopt three most prominent CSA systems (variety of crops, land levelling using laser, and zero tillage). These values were estimated as +1402, +3037 and –1577 INR ha-1, respectively, for rice-wheat cultivation system. By improving their varieties in terms of crop production, the study results showed that the farmers of the Indo-Gangetic Plains can have their net return increase to up to INR 15,712 per-hectare, per-year. Similarly, when cultivating wheat and rice with no tillage, farmers could make up to INR 6951 per-hectare, per-year, and INR 8119 per-hectare, per-year with laser-based land levelling. Given the analyses of this results, it implies that integrating individual systems together would result in an even higher yield as well as income for the farmers. In econometric terms, adoption and execution of CSA practices for crop production in the north Indian River plain would significantly influence the cost of production, which decreased, but produced an increased yield of rice and wheat.

Scherr et al. [78] reported that CSA offers to rebrand villages by providing them with embrace 'climate-smart landscapes'. This means that integrated landscape management principles that adopts the pillars of climate change must be in place prior to agricultural land allocation. The development of CSA objectives also requires strong institutional mechanism. When such systems are in place, its effects transcends to other parts of the village. Steenwerth [74] noted that while smart village residents might consider migrating to big cities, climate-smart agriculture could cause a rethink, as it gives room for entrepreneurial development in the agricultural sector, as seen in the case of youth training embarked upon in rural areas across Europe. Additionally, CSA also caters for increased demand for food due to the world's growing population. This is achieved through methods that do not jeopardize environmental health [74]. With respect to animal husbandry, some zoonotic diseases can be detected early, so that treatment plans are set underway to prevent the farmer from infection. CSA also motivates the achievement of sustainable development goals through agricultural practices that use techniques that can drive food security, improve resilience, and effectively manage emissions [70]. CSA practices are also able to curb environmental challenges related to water pollution through the use of agrochemicals [91]. A notable aspect where smart agriculture surpasses expectations is the possibility of using it as a tool for enterprise resource planning, through which the safety of agricultural products/foods can be monitored [46].

## **4. Discussion**

## *Revisiting the Research Question*

How will smart- and climate-smart agricultural research give account of, and conceptualize transformation and change in the smart village context?

In responding to the modified research question above, it is important to draw important ideas from the definitions of smart- and CSA. Albeit, CSA bears all characteristics of smart agriculture, with a step further in lowering GHG emissions. Consequently, accounts of conceptualizing transformation and change in smart village context might tend towards the adoption of key aspects of climate-smart agriculture (see Figure 1), which are somewhat multi-disciplinary in nature [92]. What this implies is that for smart villages to reach desired level in terms of development through research and policy frameworks, ideas of climate-smartness must be fully embedded across the facets of smart village agenda. According to Katara et al. [93], continuous adoption of new technologies is the first way to conceptualize transformation of smart villages. Since technology is bound to continually change, it becomes easy to bring evolving and smarter changes to smart village progress. This means that rural population must fully embrace ICT, especially since smart village idea is based on the fact that technology is adopted to hasten the growth of sustainable development [93]. Secondly, efficiency and productivity are not completely new words in smart village research. Nevertheless, it might be useful for smart village policy analysts to learn from prevention of losses for which CSA is known [94]. Another important aspect through which transformation of smart villages can be conceptualized is through capacity building of rural dwellers. As in the case of climate-smart agriculture, building capacities would bring about self-sufficiency for persons within these communities, thus reducing urban migration [15,19]. This is part of the current efforts within the different smart village initiatives. Of all the initiatives that smart village research can draw from climate-smart agricultural practices, the idea of seeking and promoting "new knowledge" [95] might be technically referred to as the most significant. Given that the world has now embraced a knowledge-based economy for which smart village development has to be a part.

On the basis of towing a part of steady development in its processes, future smart village developmental projects need to adopt successful projects of the past as a yardstick for planning. For instance, tremendous success was recorded by the IEEE smart village initiative; the EU smart village-drive, as well as the CCACFS projects, to mention a few. By adopting the recipe for success within these projects, more smart-village projects would be actualized in many parts of the world. Furthermore, it is noteworthy to state that existing smart village projects also have unique challenges. Notable amongst the challenges faced by smart villages within the IEEE project is the issue of maintenance and repairs [96]. Although as part of the project framework, two individuals are often selected and trained within the villages to fix damaged smart inputs, when demands for these inputs become high, the number of technicians might no longer be sufficient to cater for repair and maintenance needs. This is one aspect where smart village development must learn from climate-smart ideas where capacity development is well-planned and readily available.

Another aspect for which smart village development can gain from climate-smart agriculture is in its sustainability approach. While CSA strives for the cheapest routes to progress in agriculture, smart village development mostly depends on donations and funding, which slows down the pace of making progress and achieving sustained growth. As a result, for any smart villages project to achieve lasting success, such a project would have to plan self-funding strategies [10], where inputs within the village is used to generate income that would fund new projects for growth, rather than unduly wait for funding before progress is made. In building its growth, smart village planners might need to prioritize new knowledge and link it to new technology for early warning measures against potential environmental disasters. Furthermore, proponents also need to ensure that pillars of climate change are largely considered in building infrastructures [13]. This is because the impact of climate change might continue to be felt for a long time.

While ideas drawn from smart and climate-smart agriculture might indeed be useful for smart village development, Hargreaves et al. [97] explained that specific policies grounded in the values of rural areas are needed to help them transform into smart villages. This transformation must, therefore, bring effective utilization and management of resources within smart villages. The idea of transformation within the context of smart villages mostly draws attention to digital transformation, which is very important [98]. Another result of technological change is the social changes it brings [99,100].

Given the forgone discussion on how smart village development can be spurred from ideas borrowed from smart- and climate-smart agriculture, we argued that the development of a smart village has to be a gradual process. This is because the development must systematically and strategically prioritize the most important aspects, such as clean energy management and agriculture, bearing in mind the sustainability of the process.

## **5. Current Lessons & Future Research Direction**

Overall, this study revealed a number of lessons from smart-agriculture and climate-smart agriculture ideas, which, if adopted in smart villages, would achieve the following goals.


While the above lessons are specific to smart village development, there are specific shortcomings of climate-smart agriculture that must be noted [91,92], and which ought not to be adopted in smart village development.


As a result of some of the aforementioned challenges of CSA, future studies might look at the challenges posed by the adoption of "climate-smart" agriculture, prior to the full adoption of its fundamental aspects, as described in this study. This is because research by Taylor [86] pointed out certain foundation faults in the description of CSA by the World Bank group. There is also a fundamental problem in how CSA handles climate justice [87]. Another aspect opened to future research is the development of the climate-smart villages, as used in some studies [69]. While this has been achieved in some parts of the world today [13], it might be the case that smart-village research is yet to reach a maturity level as to warrant even more terms to be coined from it.

## **6. Conclusions**

The uniqueness of smart village projects around the world means that approaches towards smart village development might also differ. This study showed that smart and CSA are key areas that must be considered in developing a smart village project, and offer several lessons to proponents of smart village ideas, given how these concepts have enjoyed steady conceptualization in the research literature. Another important consideration that must be carefully explored is the tendency of developing smart villages in line with the concepts upon which smart cities are built. Having clarified in Section 1 that smart villages are not extensions of smart cities [106], it is important to understand that the challenges of rural areas differ significantly from those of cities. Hence, smart village development must come with uniquely defined plans and strategies for its development [107].

A major driving force for "smarting up" rural areas is the mass exodus of persons to the cities, as well as inferior services offered in these villages [107]. Nevertheless, an introduction of the smart village concept comes with new opportunities brought about by technology, which is currently touted as the major economic driver of the 21st century. The current study, therefore, tries to adopt the technological ideas of CSA in creating a foundational path for smart village development. To do this, the study carefully analyzes the framework of CSA and proposes that the same be adopted for developing smart villages. It is observed that certain fundamental aspects of technological innovation; productivity, new knowledge, new technology, capacity building, vulnerability reduction, increased profits, etc., are fundamental to the building of smart villages. Nevertheless, these fundamental terms cannot be embedded immediately. Rather, it must follow gradual process that gives priority to the important aspects.

**Author Contributions:** Conceptualization, A.A., O.F., O.K., A.S. and K.K.; methodology, K.K., P.M., O.K., O.F. and A.S.; software, A.A., M.A. and O.F.; validation, A.S., K.K., P.M. and O.K.; formal analysis, O.F., K.K. and O.K.; investigation, A.A., O.F., M.A.; resources, O.K., K.K. and A.S.; data curation, A.A., O.F. and O.K.; writing—original draft preparation, A.A., O.F., M.A., K.K., P.M., A.S. and O.K.; writing—review and editing, A.A., O.F., K.K. and O.K.; visualization, A.A., O.F. and M.A.; supervision, A.S., P.M., K.K. and O.K.; project administration, O.F., O.K., P.M. and K.K.; funding acquisition, O.K. and K.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded in part by the project (2020/2205), Grant Agency of Excellence, University of Hradec Kralove, Faculty of Informatics and Management, Czech Republic; project at Universiti Teknologi Malaysia (UTM) under Research University Grant Vot-20H04, Malaysia Research University Network (MRUN) Vot 4L876 and the Fundamental Research Grant Scheme (FRGS) Vot5F073, supported under Ministry of Education Malaysia for the completion of the research.

**Conflicts of Interest:** Authors declare no conflict of interest.

## **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Urban Green Infrastructure Monitoring Using Remote Sensing from Integrated Visible and Thermal Infrared Cameras Mounted on a Moving Vehicle**

**Sigfredo Fuentes \*, Eden Tongson and Claudia Gonzalez Viejo**

Digital Agriculture, Food and Wine Sciences Group, Faculty of Veterinary and Agricultural Sciences, School of Agriculture and Food, Parkville, VIC 3010, Australia; eden.tongson@unimelb.edu.au (E.T.); cgonzalez2@unimelb.edu.au (C.G.V.)

**\*** Correspondence: sfuentes@unimelb.edu.au

**Abstract:** Climate change forecasts higher temperatures in urban environments worsening the urban heat island effect (UHI). Green infrastructure (GI) in cities could reduce the UHI by regulating and reducing ambient temperatures. Forest cities (i.e., Melbourne, Australia) aimed for large-scale planting of trees to adapt to climate change in the next decade. Therefore, monitoring cities' green infrastructure requires close assessment of growth and water status at the tree-by-tree resolution for its proper maintenance and needs to be automated and efficient. This project proposed a novel monitoring system using an integrated visible and infrared thermal camera mounted on top of moving vehicles. Automated computer vision algorithms were used to analyze data gathered at an Elm trees avenue in the city of Melbourne, Australia (*n* = 172 trees) to obtain tree growth in the form of effective leaf area index (*LAIe*) and tree water stress index (TWSI), among other parameters. Results showed the tree-by-tree variation of trees monitored (5.04 km) between 2016–2017. The growth and water stress parameters obtained were mapped using customized codes and corresponded with weather trends and urban management. The proposed urban tree monitoring system could be a useful tool for city planning and GI monitoring, which can graphically show the diurnal, spatial, and temporal patterns of change of *LAIe* and TWSI to monitor the effects of climate change on the GI of cities.

**Keywords:** urban tree management; tree monitoring; computer vision; tree water stress index; leaf area index

## **1. Introduction**

Green infrastructure (GI) has become a priority in most cities worldwide and has been recognized as an essential element in urban planning and development. The urban green infrastructure, which includes natural vegetation, parks, street trees, green roofs, and small gardens, provides various benefits to the environment, community, and the economy. The GI of a city contributes valuable benefits such as regulation and reduction of temperature during heatwaves [1] through plant transpiration [2], while green roofing decreases albedo [3]. GI improves air quality [2] and reduces flood risk [4,5] and stormwater pollution [6], among other environmental benefits. Beyond improving the ecosystem, GI has been linked to improving people's physical and mental health [1]. Within the major challenges in urban cities are extreme heat and the urban heat island effect (UHI). UHI mainly occurs within cities with a higher proportion of concrete in relation to their green infrastructure (GI) [7]. In these cities, the ambient temperature increase can be multiplied by a 1.4–15 factor depending on circumstances within and surrounding a particular city environment [8,9]. UHI may worsen in the future, corresponding to the predicted climate change.

The maintenance of GI in cities may pose a challenge for city councils due to the number and complexity of tree and plant species, especially in cities classified as forest

**Citation:** Fuentes, S.; Tongson, E.; Gonzalez Viejo, C. Urban Green Infrastructure Monitoring Using Remote Sensing from Integrated Visible and Thermal Infrared Cameras Mounted on a Moving Vehicle. *Sensors* **2021**, *21*, 295. https://doi.org/10.3390/s21010295

Received: 26 October 2020 Accepted: 30 December 2020 Published: 4 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

cities, such as Melbourne, Australia [1,2]. The City of Melbourne established a goal as part of the Urban Forestry Strategy to plant 3000 trees per year as one of the primary strategies for the climate adaptation program. The city council aimed to increase canopy to 140,000 trees in 2040, twice the coverage of the existing number of trees that it manages at present [3]. Trained arborists maintain public trees through routine inspection and assessment, which can also be requested through public reporting. Trees situated in heavily accessed areas such as parks and boulevards are inspected annually, while other locations are inspected at least once in 2 years. The city eliminates about 800 tree stands due to various reasons, including threats to safety. A considerable population of trees managed by the City of Melbourne are a century old and may pose higher risk of decline and, consequently, safety risk.

The manual inspection of urban trees by arborists is quite time demanding and inefficient. With the large-scale greening plans in urban cities such as Melbourne, it is highly impractical and nearly impossible to monitor GI to achieve high temporal and spatial resolution through the current manual practice. Other methods use wireless sensor networks, including the monitoring of soil moisture, temperature, light, humidity, and pressure [10], using internet of things (IoT) systems [4], real-time controls [5], and plant/tree-based sensors, such as sap flow probes [6,11]. These methods produce accurate and high temporal data resolution. However, one major disadvantage is that this method can only be applied to a few representative plants, as installing sensors to numerous trees is costly. Furthermore, the sensor networks require frequent monitoring and high maintenance, requiring specialized personnel with specific technical skills.

More spatially representative approaches for GI monitoring of cities are based on satellite remote sensing on relevant vegetation indices (VIs) related to growth and water status [12]. This has been applied in China for 70 major cities [13], and in Sweden [14] using Sentinel-2 and Landsat-8 satellites [13], as well as in Croatia using World View 1, 2, and 3 with high resolution visible and multispectral bands [15]. The use of satellite imagery can monitor large areas from a single image or stitched up images incorporating several square kilometers, which is a major advantage. Also, information can be readily available, and sources can be either free or low-cost (i.e., Landsat and Sentinel satellites). However, disadvantages can include low resolution of information per pixel, reaching 0.5 m for panchromatic imagery, and between 2 m to 30 m per pixel for multispectral imagery. Higher spatial resolution imagery may be expensive, such as those from the World View satellites. Furthermore, satellite revisit time to the same spot (i.e., in cities) may be between 10–15 days, and data quality depends on how clear the skies are.

To address the problem of low temporal and spatial resolution, airborne, and unmanned aerial vehicles (UAV) have been implemented to monitor GI in cities [16–19]. However, the use of airborne remote sensing comes with a cost, requiring a pilot and skilled personnel to operate the instrumentation, process the information, and deliver interpreted information to relevant city council personnel for GI management and decision making. Some services such as Nearmap (Nearmap, Barangaroo, NSW, Australia) offer high-resolution visible images with a high temporal resolution for major cities and coastlines [20,21]. However, the application of visible images is mainly for monitoring of growth parameters for trees [22]. The recent popularity of UAV has also expanded its application in remote sensing, with the accommodation of various camera and sensor payloads aside from visible such as multispectral camera and LIDAR [23–26]. With UAV, the main challenge is the implementation in countries with strict civil aviation regulations, such as many cities in Europe, the United States, and Australia, among others. In Australia, for example, the Civil Aviation Safety Authority (CASA) has a very strict regulation to fly drones within 30 m in proximity to people, making flights in heavily populated cities, such as Melbourne (Victoria, Australia) virtually impossible [27].

This paper proposes a novel GI monitoring approach based on prototype integrated visible and infrared thermal cameras to automatically obtain different VIs based on growth and tree water status parameters on a tree-by-tree scale. The integrated cameras are

mounted on top of moving vehicles, circulated through one of the most important and historical Elm tree avenues in Melbourne (Royal Parade), Australia. The integrated system, which is composed of low-cost instrumentation, could be mounted on top of public transport vehicles, such as buses and trams, city council vehicles, and rubbish trucks. Public transport vehicles allow for incursion to potentially every street at multiple times in a day, offering a high temporal and spatial resolution of trees monitoring. The proposed system enables automated data acquisition, analysis, and mapping and does not require specialized personnel. It can also offer diurnal monitoring of trees to assess in real-time the effects of weather anomalies, such as heatwaves, floods, and heavy winds, among others, and the detection of pest and disease incidence. The novel technology (integrated cameras) and application could potentially be an accurate, cost-effective, and user-friendly tool for city councils, and they could base management strategies with high reliability on the system proposed, such as tree lopping, detection of encroachment of tree branches on power lines, and deterioration of old trees.

Only a few studies are based on the implementation of cameras for upward-looking imagery and analysis without any automation. These have been restricted to the 3D modeling of trees using handheld cameras and point cloud analysis [28] and the analysis of trees' thermal characteristics [29]. Other applications using different technologies, such as hyperspectral cameras [30] and low-cost electronic noses (e-noses) [31], using the methodology proposed can be implemented to obtain more information from trees and their environment, such as diagnosis of vegetation health.

This study has been based on the integration of previously developed technology from our research group for the automated analysis of visible and infrared thermal imagery for different crops such as eucalyptus trees [22], grapevines [32–35], kiwi plants [36], apple trees [37], cherry trees [38–40], and cocoa plants [41], among others.

## **2. Materials and Methods**

Data acquisition was performed in Melbourne, Australia, mainly based on an integrated visible and infrared thermal camera developed and processed using customized computer vision algorithms.

## *2.1. Urban Site and Tree Material Description*

The monitoring site (Figure 1) was located along the iconic Royal Parade avenue in the city of Melbourne, Australia, that starts at Grattan Street (−37◦48 02.27 S; 144◦57 26.27 E; 33 m.a.s.l.) finishing on Park Street (−37◦46 41.45 S; 144◦57 36.56 E; 46 m.a.s.l.), and vice versa. The trees are planted along a nature strip separating the main road and the access road in both directions. The main roads (North and South bound) are divided by a median strip containing the Route 19 tram lane. Each way corresponds to 2.52 km, with a total distance of 5.04 km both ways. There are 172 deciduous trees considered in the monitoring, composed of different Elm species (*Ulmus* spp.) planted in 1900 and 1997 [42]. The trees are irrigated using sub-surface irrigation, and tree lopping management is performed regularly by the Melbourne city council.

## *2.2. Climate and Weather Information Description*

The climate in Melbourne is classified as subtropical oceanic with mild winters and pleasant to hot summers. Windy conditions are common, and weather changes can occur within the same day. The average temperatures between November and January are between 22 and 26 ◦C, with minimum temperatures between 11 and 14 ◦C. The yearly average precipitation is 670 mm, with an even distribution throughout the year of around 50 mm per month. Sunshine hours are higher between September and March (between 6–9 h). Specific weather data available for the trial site and the monitoring period were acquired from the Bureau of Meteorology, measured from a meteorological station located in Melbourne Olympic Park (Number: 086338) at 3.4 km from Royal Parade. The weather information extracted from this station was: maximum daily temperature (◦C), rain (mm),

and solar radiation (MJ m<sup>−</sup>2). Monitoring was performed from November 2016 (late spring) to January 2017 (summer).

**Figure 1.** Location monitored using a moving vehicle from (**A**) start of Royal Parade from Grattan Street to (**B**) Park Street (2.52 km), and vice versa. There were 172 trees monitored, consisting of different species of Elm trees (*Ulmus* spp.).

## *2.3. Integrated Visible and Thermal Infrared Camera System*

The integrated camera system (Figure 2) consisted of a visible RGB video camera and a thermal infrared camera FLIR AX8™ (FLIR Systems, Wilsonville, OR, USA) with a resolution of 90 × 60 pixels, connected to a web-based system that can simultaneously capture and store the videos and infrared thermal images (IRTIs) to be further downloaded for analysis or transmitted to cloud storage and processing system. The thermal camera had a spectral range of 7.5–13 μm, an accuracy of ±2 ◦C, and an emissivity of 0.985. The IRTI capture rate was every second. The RGB video camera is connected to a Raspberry Pi Camera Module V2.1 (Raspberry Pi Foundation, Cambridge, UK; Figure 2A), board, and memory card. This device has an 8-megapixel sensor with a resolution of 640 × 360 pixels, 4:3 aspect ratio, and 30 frames per second (fps). Videos were recorded within the unit in H.264 video compression format and automatically converted into Motion Pictures Expert Group-4 (.mp4) files. The camera was fitted with a 3-axis gimbal to minimize movements when acquiring the data (Figure 2A); an integrated temperature, relative humidity, and solar radiation sensors within a 3D printed Stevenson screen (Figure 2A); and a magnetic GPS tracker (Figure 2B). The integrated camera was mounted on top of a car (Figure 2B) with a height between the camera and the tree canopies of approximately 5 m.

**Figure 2.** The integrated camera system. (**A**) The system is composed of weather and shock resistant case (1) that holds the Raspberry Pi boards and battery; thermal infrared camera FLIR AX8™ (2); the visible Red, Green, and Blue (RGB) Raspberry Pi Camera Module V2.1 (3); power mount receptacle (4) to charge the internal battery; 3-axis gimbal (5) to provide stability to the camera; integrated temperature, relative humidity, and solar radiation sensors (Stevenson screen, 6). (**B**) Example of the mounting procedure of the integrated camera on top of a vehicle. The camera was also integrated with a magnetic GPS tracker (7).

The radiometric data from the thermal infrared camera were obtained every second while traveling through the Royal Parade. They were recorded as in comma-separated values (.csv) format files and the visible RGB images in Joint Photographic Experts Group (.jpg). Both sets of data were obtained using the Sense Batch software (SENSE Software, Warszawa, Mazowsze, Poland). The data were analyzed using customized codes developed and updated using Matlab® R2020b (Mathworks Inc., Natick, MA, USA).

The camera's integrated sensors consisted of an AM2302 (wired DHT22) temperaturehumidity sensor (Guangzhou Aosong Electronics Co., Ltd., Guangzhou, China). This sensor can obtain new data from it once every 2 s (0.5 Hz), which is accurate for 0–100% humidity readings with 2–5% accuracy and −40 to 80 ◦C temperature readings with ±0.5 ◦C accuracy. The SP-510-SS upward-Looking Thermopile Pyranometer (Apogee Instruments, Inc., Logan, UT, USA) has a sensitivity of 0.05 mV per W m<sup>−</sup>2, with a measurement range between 0 to 2000 W m−<sup>2</sup> (net shortwave irradiance) and repeatability of <1%. The detector response time is 0.5 s with a field of view of 180◦ and spectral range of 385–2105 nm, directional (Cosine) response less than 30 W m−<sup>2</sup> at 80◦ solar zenith, temperature response: <5% from −15 to 45 ◦C at the operating environment: −50 to 80 ◦C, and 0 to 100% relative humidity.

#### *2.4. Image Pre-Processing and Computer Vision Algorithms*

Every frame corresponding to a canopy from the visible (RGB) video and infrared thermal images were analyzed using the computer vision algorithms described in Sections 2.4.1 and 2.4.2, respectively. Figure 3 shows an example of a visible (RGB) frame and corresponding infrared thermal image from an Elm tree canopy along the Royal Parade.

**Figure 3.** Example of a visible (RGB) image (**A**) and the corresponding infrared thermal image (**B**) from an Elm tree taken using the integrated camera on top of a moving vehicle.

The pre-processing of the RGB images consisted of the binarization (Figure 4A) using the blue channel from the RGB images by selecting the lowest part of the histogram curve (valley) detected automatically between the pixels corresponding to the canopy material (first peak) and the background or sky (second peak) (Figure 4B). After binarization, each image was automatically subdivided into a 5 × 5 sub-images to perform gap analysis. A large gap (lg) per sub-image was considered when there was over 75% of sky. Total pixels (tp) corresponded to a fixed value related to the resolution of the camera used. This pre-analysis has been described in detail in Fuentes et al. [22].

The pre-processing of the thermal images was performed in batch after each measurement campaign using the SENSE Batch software (Sense Software, Warszawa, Mazowsze, Poland), which extracts radiometric data per pixel in a comma-separated file (.csv) in the form of a matrix processed in Matlab (Figure 4C). Leaf material was selected by simple automatic elimination of temperatures below 0 + ◦C since this separates the sky from the canopy material (Figure 4D). From the segmented image, the canopy temperature was automatically extracted (T*canopy*) as entry parameter for the TWSI and Ig calculation (Equations (7) and (8)).

#### 2.4.1. Canopy Architecture and Growth Parameters

Videos from the visible camera were processed automatically using a customized code written in Matlab® R2020b to analyze frames following a computational process proposed by Fuentes et al. (2008) [22].

Canopy architecture parameters were obtained using the following algorithms considering the fractions of foliage projective cover (*ff*), crown cover (*fc*), and crown porosity (*Φ*), which were calculated using the following computational algorithms proposed by Fuentes et al. (2008) [22]:

$$f\_f = 1 - \frac{t\mathcal{g}}{tp} \tag{1}$$

$$f\_{\mathfrak{c}} = 1 - \frac{l\_{\mathfrak{g}}}{tp} \tag{2}$$

$$\phi = 1 - \frac{f\_f}{f\_c} \,\tag{3}$$

where *lg* = large gap pixels, *tg* = total pixels in all gaps, and *tp* = total gap pixels.

**Figure 4.** Example of the automated pre-processing of visible (RGB) images transformed to binary images (**A**) using the blue channel as filter (**B**) and the corresponding infrared thermal image (**C**) filtered to create a mask (**D**) to account for leaf material.

*LAI* (adimensional) is calculated from Beer's Law, defined as the total one-sided area of leaf tissue per unit 3 ground surface area [43]. Hence, the *LAI* values describe m2 of leaf area per m2 of soil.

$$LAI = -f\_c \frac{\ln \phi}{k} \tag{4}$$

where *k* = coefficient of light extinction (*k* = 0.5), which is applicable for tall trees [22], and the clumping index at the zenith, Ω(0), was calculated as follows:

$$
\Omega(0) = \frac{(1 - \phi)\ln\left(1 - f\_f\right)}{\ln(\phi)/f\_f}.\tag{5}
$$

The clumping index is a correction factor in obtaining effective *LAI* (*LAIe*), also adimensional, which is the product of:

$$LAI\_{\mathfrak{c}} = LAIx\Omega(0). \tag{6}$$

Equation (5) describes the non-random distribution of canopy elements. If Ω(0) = 1 means that the canopy displays random dispersion, then for Ω(0)> or <1, the canopy is defined as clumped.

## 2.4.2. Infrared Thermal Image Analysis

A tree water stress index (TWSI) was derived from the common crop water stress index (CWSI) [32] used in agriculture, which is a normalized value (0–1) and, therefore adimensional, and it was calculated using the following equation after determining *Tdry* and *Twet* [44]:

$$\text{TWSI} = \frac{T\_{canopy} - T\_{wet}}{T\_{dry} - T\_{wet}} \tag{7}$$

where *Tcanopy* is the actual canopy temperature extracted from the thermal image at determined positions, and *Tdry* and *Twet* are the reference temperatures (in ◦C) obtained using the statistical temperature distribution discrimination described in published research [39].

An infrared index (*Ig*), which is adimensional and proportional to leaf conductance and water vapor transfer (*gs*), can be obtained using the relationship as follows [45]:

$$I\_{\rm g} = \frac{T\_{\rm campy} - T\_{\rm uvt}}{T\_{\rm dry} - T\_{\rm uvet}} = \ g\_s \left( r\_{\rm av} + \left( \frac{s}{\gamma} \right)\_{r\_{\rm HR}} \right) \tag{8}$$

where *raw* = boundary layer resistance to water vapor, *γ* = psychrometric constant, and *s* = slope of the curve relating saturation vapor pressure to temperature [45,46].

For automated analysis, the leaf energy balance approached was implemented using integrated sensors within the camera described in Figure 2A as [32]:

$$Tdry - Ta = \frac{r\_{HR}R\_{wi}}{\rho c\_p} \tag{9}$$

where *Ta* is the air temperature measured at the same positions and time as infrared thermography acquisition, *rRH* = the parallel resistance to heat and radiative transfer, *Rni* is the net isothermal radiation (the net radiation that would be received by an equivalent surface at air temperature), *ρ* is the density of air, and *cp* is the specific heat capacity of air. This formula uses the concept of isothermal radiation and assumes a dry surface with the same aerodynamic and radiative properties, in which the sensible heat loss will equal the net radiation absorbed [47].

$$Twet-Ta = \frac{r\_{HR}r\_{aW}\gamma R\_{ni}}{\rho c\_p[\gamma(r\_{aW}) + sr\_{HR}]} - \frac{r\_{HR}\delta e}{\gamma(r\_{aW}) + sr\_{HR}}\tag{10}$$

The thresholds *Twet* and *Tdry* are references that can be leaves painted with water (*Twet*) and use petroleum jelly (*Tdry*) to obtain through infrared thermography the maximum and minimum temperatures to be found within a specific canopy at the time of measurements [32]. The leaf energy balance approach allows the implementation of an automated procedure to obtain these thresholds using the sensors incorporated in the integrated camera proposed (Figure 2).

#### *2.5. Survey, Automated Detection of Trees Location, Data Extraction, and Mapping*

Acquisition of images was performed on four dates: twice in November 2016 (17 and 19 November), followed by 19 December 2016 and 16 January 2017. The image surveys were all performed at 1–2 pm during maximum atmospheric demand (maximum vapor pressure deficit), a common practice in agriculture to assess plant water status for irrigation assessment requirements.

For the 172 Elm trees monitored in this study, the GPS location was extracted from Google Earth Pro (Googleplex, Mountain View, CA, USA). The tree positions were used as anchors to automatically extract information from procedures previously explained for canopy architecture and infrared thermal-based parameters. The automated extraction consisted of identifying the nearest coordinates registered in the integrated camera to the anchored GPS for specific trees.

Once the data were extracted, they were mapped using a customized code written in Matlab® R2020b to produce: (i) geo-located icons (circles) with relative sizes to denote changes in growth (*LAIe*), and (ii) geo-located circles with a different color to represent different TWSI values. The process can be used to map any parameter extracted using Equations (1)–(10).

### **3. Results**

### *3.1. Weather Data within the Period of Measurement and Calculated Parameters*

Figure 5 shows the weather information acquired from the closest meteorological station from the trial site. The first two dates of measurement (A: 17 November 2016 and B: 29 November 2016) had maximum rain events of 12.6 and 17 mm of rain in the previous week, and maximum temperatures of 31.4 and 20 ◦C, respectively. These dates had high solar radiation (28.8 MJ m−<sup>2</sup> and 31 MJ m−2, respectively). In the last two measurement surveys (19 December 2016 and 16 January 2017), there were no or minimal rain events (0 mm and 4 mm, respectively) within two weeks preceding the measurements. Both dates had high maximum temperatures and solar radiation values (30.02 and 32.7 ◦C and 29.3 and 30.5 MJ m<sup>−</sup>2, respectively).

**Figure 5.** Meteorological data showing daily maximum temperature (◦C) solar radiation (MJ m−2) and rain (mm) in the Royal Parade for four different dates studied: (**A**) 17 November 2016, (**B**) 29 November 2016, (**C**) 19 December 2016, and (**D**) 16 January 2017.

Table 1 shows the main canopy and tree water status parameters for all the measurement survey days. There was considerable variation in *LAI* and *LAIe* from a minimum of 0.61 and 0.41, found the last date of measurement, to maximum values of 5.98 for *LAI* in the first date of measurement and 4.97 *LAIe* for the second date of measurement. Furthermore, the lowest Tc values, TD, and TWSI corresponded to the second measurement date.

**Table 1.** Growth and water stress parameters obtained using the proposed urban tree monitoring system, measured at Royal Parade in Melbourne, Australia, for four measurement surveys between 2016 and 2017. Parameters are presented with maximum, minimum, means, and standard deviation values (SD) for leaf area index (*LAI*, adimensional), effective *LAI* (*LAIe*, adimensional), canopy temperature of trees (Tc, ◦C), temperature depression (TD, ◦C), thermal infrared index (*Ig*, adimensional), and tree water stress index (TWSI, adimensional).


## *3.2. Comparative Analysis of Main Extracted Parameters from Trees*

Figure 6A compares growth parameters (*LAIe*) for the 172 trees monitored with the TWSI for the different measurement dates. The trends followed apparent curvilinear relationships with the last two dates (19 December 2016 and 16 January 2017) with lower *LAIe* and higher TWSI than the earliest dates (17 November 2016 and 29 November 2016). Figure 6B shows the comparison between the Ig and TD parameters related to stomatal conductance from trees. There was contrasting behavior of these parameters for the first two dates with lower Ig and higher TD for the first and flat distribution of the whole range of Ig values with low TD close to the 0 values. On the contrary, the last two dates had similar behavior with low Ig and TD values ranging from −3 to around 3 ◦C.

**Figure 6.** Comparison between effective leaf area index (*LAIe*, dimensionless) and tree water stress index (TWSI) (**A**), and between the infrared thermal index (*Ig*) and temperature depression (TD, ◦C) (**B**) for 172 elm trees monitored along the Royal Parade in Melbourne, Australia for four different dates between 2–16 and 2017, using the proposed urban tree monitoring system.

## *3.3. Main Growth and Tree Water Stress Parameters Map*

Figures 7 and 8 show the proposed urban tree monitoring system's main outputs, displaying the main parameters extracted per tree along Royal Parade in four measurement dates. Figure 7 shows the *LAIe* for different trees with the relative size of circles corresponding to trees changing according to growth differences between dates. Figure 8 shows changes in color of circles representing the trees relative to the TWSI for different dates.

**Figure 7.** Mapping of effective leaf area index (*LAIe*) along the Royal Parade (5.04 Km) of 172 trees using the proposed urban tree monitoring system for four different dates: (**A**) 17 November 2016, (**B**) 29 November 2016, (**C**) 19 December 2016, and (**D**) 16 January 2017. Different colors and relative circle sizes correspond to the *LAIe* scale.

**Figure 8.** Mapping of tree water stress index (TWSI) along the Royal Parade (5.04 Km) of 172 trees using the proposed urban tree monitoring system for four different dates: (**A**) 17 November 2016, (**B**) 29 November 2016, (**C**) 19 December 2016, and (**D**) 16 January 2017. Different colors correspond to the TWSI scale.

## **4. Discussion**

The proposed urban tree monitoring system that uses an integrated camera on moving vehicles can automatically provide information on trees' growth and water status changes, which can serve as a powerful decision-making tool for city councils for tree management (i.e., supply water requirement at appropriate times and tree lopping for power lines encroachment and public safety management). The reliability of the system is based on the growth and canopy architecture parameters and algorithms used, which have been successfully implemented for other trees such as eucalyptus [22], and tree crops such as cherry trees [38,40], apple trees [37], and grapevines [33,48–50]. Tree water stress algorithms have been used to describe the water status of many trees and crops [32,51–53]. Furthermore, most of the tree canopies were visible in the field of view of canopies for both visible and infrared thermal images (Figures 3 and 4), making the analysis representative of the whole tree. Furthermore, since images are upward-looking, the monitored parts of the trees were the under canopy and were shaded, which has been regarded as the most consistent and representative part to monitor using infrared thermal imagery [46,54].

The sensitivity of the growth and physiological parameters obtained and their variations are specifically shown in Figure 6 and compared between the trees measured (individually) and temporally (within dates). The variation is sensible in response to weather conditions and changes related to atmospheric demand (temperature) and water availability (rain). These trends and their sensitivity are further supported by the mapping of the processed data in the form of *LAI* (Figure 7) and TWSI (Figure 8). The parameters obtained are in accordance with weather information acquired within the measurement dates. The first two dates (17 November 2016 and 29 November 2016) corresponded to milder weather, with cooler weather for the second date with a maximum temperature of 20 ◦C, followed by rain events. For the second date, lower atmospheric demands produced a flat response for TWSI and TD. In the case of TD, lower TD values, close to 0 ◦C, are related to low stomata opening and transpiration (Figure 6B). However, they were not associated with higher TWSI (Figure 6A). The rest of the dates (19 December 2016 and 16 January 2017) have more significant increases in TWSI with higher atmospheric demands (evapotranspiration), as shown by higher maximum temperatures and solar radiation. The highest and more significant determination coefficient (Figure 6B) between Ig and TD was found for the dates with higher atmospheric demand (first, third, and fourth dates), which was expected since these parameters are related to stomata aperture [32,54].

Another advantage of the proposed system is that it allows the automatic mapping of data obtained from surveys on a tree-per-tree scale (Figures 7 and 8). For growth parameters, such as *LAIe* (Figure 7), some of the trees with higher growth showed decreased *LAIe* from the first to the second date of measurement, which may be related to continuous tree lopping management from the council (Figure 7A,B). However, the lowest and most consistent *LAIe* values were found in the last date of measurement (January 2017), which corresponded to one of the hottest months in summer and the starting of the senescence stage for the Elm trees (Figure 7D), in comparison to the previous dates. For TWSI, the parameter trends followed water availability from rain events during the last weeks and maximum temperatures with the highest values corresponding to the warmest dates in December 2016 and January 2017 (Figure 8C,D).

The integrated cameras could be mounted on public transport of cities, such as buses and trams. The installation of the system on trams is ideal, being on rails are on a fixed route, which can offer more precise data acquisition and more reliable comparative analysis. Furthermore, at least along the Royal Parade route (Route 19), a particular tram can pass through the same spot every 80–90 min, which can acquire at least 13 data points in one day from 5 am to midnight. To access more places within the city, such as suburb streets, cameras could be installed on rubbish trucks and buses, which have more extensive access to residential areas. This layout of the trams' path is similar to many European cities since they have similar designs.

The diurnal data collected may be more relevant for infrared thermal parameters to assess tree water status changes throughout the day compared to changes in growth, which are expected to be minimal, while data related to changes in leaf or branch angle due to water stress after sunset could be relevant to assess night-time water loss by trees, as this phenomenon is relevant to other tree species and crops [55–59]. Continuous daily data of water stress may also offer insights of tree behavior within heat waves [60], pest and disease interactions [61,62], windy days, and mortality estimates [63]. The volume of data that can be gathered through the proposed system allows the implementation of machine learning modeling and artificial intelligence to promptly detect problems for management and mitigation, avoiding damage to infrastructure and the public due to unpredicted fallen trees or big branches.

The system has been proposed, and data analysis can be deployed as a user-friendly digital platform producing maps with tree water status and growth maps depicted in Figures 7 and 8. Users can click on any individual tree and obtain numerical and other management information as it is already set up for planting date and basic information of trees by the Melbourne city council [3]. Furthermore, since the conception of the integrated camera idea, on which this paper was based, FLIR has released an integrated visible 4K video camera and a high resolution infrared thermal imaging: FLIR Duo Pro (FLIR Systems, Wilsonville, OR, USA). This camera is intended to be mounted as a payload for UAV vehicles. It can also be used to acquire data to obtain the analysis proposed in this paper mounted on vehicles as per Figure 2B. The downside will be the costs of using these cameras if many vehicles are required for this purpose. It is thought that the higher resolution from the FLIR camera will not impact with statistical significance results obtained with the low-cost camera system presented in this paper. The latter is supported by previous research that has compared different resolutions of visible and thermal infrared cameras for growth and water status assessment on trees with no significant differences for the parameters studied [39,64], which can be explained by the short height between the camera and the canopies included for these type of studies, which is between 3–5 m.

This study was based on an extensive avenue in which there were a predominant tree species. Hence, further studies should be conducted for different tree species to account for the variety that exists in a normal urban green infrastructure environment. Even though the algorithms used in this study have been proven to be robust for other horticultural tree species, specific calibrations should be made to consider different canopy architectures and sensitivity/tolerance to different water stress levels.

#### **5. Conclusions**

The urban green infrastructure could be automatically monitored using a low-cost integrated camera system mounted on top of moving vehicles. Specifically, the main advantages of the system described in this paper compared to similar studies to monitor the green infrastructure in urban environments are: (i) low-cost instrumentation required to integrate visible and infrared thermal cameras; (ii) the system can be mounted on public transport such as buses, trams, and city council vehicles with the extra advantage when considering garbage trucks since they can access every street of a city if extensive monitoring is required; (iii) it could provide high spatial and temporal data resolution, which is related to the frequency of public transport through the same trees; (iv) algorithms implemented are robust and have been successfully tested on a wide variety of horticultural trees; (v) the system does not require special permits or trained pilots, such as the case of UAVs, and they also do not have restrictions due to privacy issues since they monitor urban infrastructure in an upward-looking fashion above the pedestrian level. These operational, cost-effectiveness, accuracy, and privacy-related advantages of the system proposed can be compared to those of manual measurements of green infrastructure, using sensors and IoT on sentinel trees, remote sensing using satellites, UAVs, or the airborne instrumentation (Nearmap) discussed in this paper. Furthermore, the high volume of data collected (spatial and temporal) using the system proposed in this paper could allow the implementation

of machine learning algorithms and artificial intelligence (AI) to obtain further vegetation indices of trees to manage the cities' green infrastructure efficiently, to maximize resources, and to minimize detrimental effects of climate change and risk to infrastructure and people.

**Author Contributions:** Conceptualization, S.F. and E.T.; methodology, S.F.; software, S.F.; validation, S.F., E.T. and C.G.V.; formal analysis, S.F., E.T. and C.G.V.; investigation, S.F.; resources, S.F. and E.T.; data curation, S.F., E.T. and C.G.V.; writing—original draft preparation, S.F.; writing—review and editing, E.T. and C.G.V.; visualization, S.F.; funding acquisition, S.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially funded by the Melbourne Network Society, belonging to The University of Melbourne.

**Institutional Review Board Statement:** Not Applicable.

**Informed Consent Statement:** Not Applicable.

**Data Availability Statement:** Data and intellectual property belong to The University of Melbourne; any sharing needs to be evaluated and approved by the University.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


## *Article* **Machine-Learning Classification of a Number of Contaminant Sources in an Urban Water Network**

**Ivana Luˇcin 1,2,\*, Luka Grbˇci´c 1,2, Zoran Carija ˇ 1,2 and Lado Kranjˇcevi´c 1,2**


**Abstract:** In the case of a contamination event in water distribution networks, several studies have considered different methods to determine contamination scenario information. It would be greatly beneficial to know the exact number of contaminant injection locations since some methods can only be applied in the case of a single injection location and others have greater efficiency. In this work, the Neural Network and Random Forest classifying algorithms are used to predict the number of contaminant injection locations. The prediction model is trained with data obtained from simulated contamination event scenarios with random injection starting time, duration, concentration value, and the number of injection locations which varies from 1 to 4. Classification is made to determine if single or multiple injection locations occurred, and to predict the exact number of injection locations. Data was obtained for two different benchmark networks, medium-sized network Net3 and largesized Richmond network. Additionally, an investigation of sensor layouts, demand uncertainty, and fuzzy sensors on model accuracy is conducted. The proposed approach shows excellent accuracy in predicting if single or multiple contaminant injections in a water supply network occurred and good accuracy for the exact number of injection locations.

**Keywords:** water distribution networks; water network contamination; machine learning; random forest; neural network

## **1. Introduction**

Contamination in water distribution networks can occur due to deliberate or unintentional intrusions and it is of extreme importance to determine the contamination event parameters so it can be detected which parts of water distribution networks have been exposed to the contaminant and needed measures can be conducted. This is considered to be an inverse problem since injection location, injection starting time, injection duration, and contaminant chemical concentration value needs to be predicted based on sensor measurements. Numerical simulations are used to determine these parameters, but model limitations need to be taken into consideration. EPANET [1] is the most commonly used software for water distribution network simulations and uses an advective approach which cannot efficiently analyze contaminant dispersion in the networks. Piazza et al. [2] conducted experiments where it was shown that dispersive and diffusive processes must be incorporated in the transport model for less turbulent fluid flows to achieve more accurate results than the pure advection model. Also, EPANET assumes complete mixing in all network junctions, which can be valid only in the case of a single outlet or if there is considerable distance between two junctions. Therefore, EPANET extension EPANET-BAM [3] was proposed which uses experimentally calibrated mixing model parameter to more accurately model mixing in network junctions. A number of studies investigated mixing behavior for different conditions, both experimentally and numerically, to further enhance these simpler 1D numerical models [4–9].

**Citation:** Luˇcin, I.; Grbˇci´c, L.; Carija, ˇ Z.; Kranjˇcevi´c, L. Machine-Learning Classification of Number of Contaminant Sources in an Urban Water Network. *Sensors* **2021**, *21*, 245. https://doi.org/10.3390/s21010245

Received: 26 November 2020 Accepted: 29 December 2020 Published: 1 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/ licenses/by/4.0/).

Huang and McBean [10] investigated a data mining approach for identifying possible sources of intrusion where single and multiple injection scenarios were considered. In the case of multiple injection scenario, the method provided a limited number of nodes with the probability of them being the true contamination source. However, in their work, it is not predicted what is the true number of injection locations. In Wang and Harrison [11] a Bayesian approach was coupled with Support Vector Regression to provide a probability distribution of water network nodes being contaminant sources. However, a single injection is assumed, and it is noted that multiple contaminant sources should be considered in future work where the likelihood evaluation needs to be adjusted. Seth et al. [12] investigated the efficiency of three different methods for source detection; Bayesian probability-based method, backtracking method (using contaminant status algorithm), and optimizationbased method where accuracy in case of multiple injection locations was investigated for two and three contamination injection locations. It was noted that the Bayesian method is designed only for a single contamination location while the contaminant status algorithm used in De Sanctis et al. [13] provides a list of possible solutions that narrow down search space for the optimization method; however, it also does not identify the possible number of injection locations. In Luˇcin et al. [14] a new search space reduction method was proposed, which can eliminate a considerable number of source nodes for both single and multiple injection locations, but with considerably greater reduction for single injection scenario. A number of different optimization approaches were considered to determine the contamination source, an overview of proposed methods can be found in Adedoja et al. [15]. Optimization approach can be easily extended to consider multiple contamination sources, as mentioned in [16–18].

If considering the optimization approach with multiple injection locations, with each additional source of contamination, the complexity of search space increases with an increase of optimization variables. Since the number of injection locations is not known, as a precaution, multiple injection locations should be allowed, since optimization can set variables to zero (which eliminates that source node and eliminates the number of injection locations), but it cannot add additional variables (injection locations) during the optimization process. In this way, in the case of a single injection location, optimization can eliminate other source nodes (all contamination parameters would be set to 0). However, this considerably increases the complexity of the considered problem since unnecessary fitness function evaluations would be conducted due to greater search space. Thus, it would be greatly beneficial to determine the number of injection locations before the optimization algorithm is employed. Also, if it is known that a single injection event occurred, a number of methods can be used more efficiently to reduce the complexity of the problem. For example, the machine learning approach provides probabilities for each network node being the true contamination scenario, which greatly reduces the number of suspect nodes and helps in quicker detection of true contamination location. However, in the case of multiple injections, different likelihood evaluation is needed which increases the complexity of the machine learning approach. Prediction of the number of contamination sources has previously been conducted for air pollution in Wade and Senocak [19], but to authors knowledge was not conducted for water distribution network contamination scenarios.

Machine learning tools have been increasingly used in contamination detection, where Random Forest has been used for groundwater source of contamination detection [20] and source detection in a river [21]. In Grbˇci´c et al. [22] Random Forest algorithm was used to predict contamination event parameters in water distribution networks and in Grbˇci´c et al. [23] new machine learning-based algorithm was proposed. A great advantage of prediction models is that they can be constructed before an accident occurs, so when a contamination event is detected prediction can be made even for large networks in a computationally efficient way. Thus, the proposed model which predicts number of injection locations can be used prior to conducting approaches that search for contamination parameters, without influencing the reaction time needed to contain the contamination event. However, in accident situations hydraulic conditions can greatly differ from those

on which model was trained, thus, a wrong prediction could be made. This can be handled with the preparation of multiple prediction models with different hydraulic conditions or by using a prediction model that achieves great accuracy with the small number of inputs so time for prediction also becomes negligible considering the benefit of search space reduction when redundant optimization parameters are not used.

In this paper, the Random Forest and Artificial Neural Network classifier are used to predict the number of contamination sources based on contamination sensor measurements in the water distribution network. Sensor measurements of contamination needed for model teaching are obtained from contamination scenarios simulated using EPANET2 with Monte Carlo generated contamination parameters. An investigation was conducted for two different sized benchmark water distribution networks with different sensor layouts, to examine the efficiency of the proposed machine learning approach. Investigation of demand uncertainty and fuzzy sensors is also estimated.

## **2. Materials and Methods**

## *2.1. Benchmark Water Supply Networks*

Prediction of the number of injection sources is conducted for two benchmark different sized networks. Investigated networks are Net3 EPANET2 example consisting of 92 nodes and Richmond network consisting of 865 nodes, obtained from The Centre for Water Systems (CWS) at the University of Exeter [24]. For the Net3 network, two different sensor layouts are investigated. In first layout four sensors were placed in network nodes 117, 143, 181, and 213 as in [25] and in second layout four sensor were placed in network nodes 115, 119, 187, and 209 as in [26]. Additionally, an investigation of the number of sensors was conducted. For the first layout, two sensors were placed in network nodes 117 and 181, and for the second layout sensors were placed in network nodes 119 and 209. For Richmond network five sensors were placed in network nodes 93, 352, 428, 600, and 672 where sensor layout was taken from [27]. Layout with three sensors placed in network nodes 93, 428, and 672 was also considered. Considered networks with sensor layouts can be seen in Figures 1 and 2.

Contamination scenarios are simulated using EPANET2 version 2.0.12. where for both networks, simulation time is 24 h with a hydraulic time step of 10 min, quality time step 5 min, pattern time step 10 min and report time step 1 h. For all conducted simulations, the EPANET2 flow paced method is used for the contaminant injection. Contamination scenario parameters are chosen randomly. The number of injection locations is chosen from 1 to 4 nodes. The starting time and duration of contamination injection are chosen from 0 to 24 h. Concentration was randomly chosen from 10 to 2000 mg/L. For contamination scenarios with multiple injection locations starting time, duration, and concentration was kept the same for every injection location.

Prior to simulating multiple injection scenario, independent simulations for each randomly chosen node as a source of contamination are conducted. If contamination is not registered for the investigated node with chosen contamination parameters, that node is eliminated as source location and only nodes for which contamination was detected in at least one sensor are kept as a source of contaminant. For example, if four source nodes are randomly chosen to be the source of contamination, but only two source nodes influence sensor detection of contaminant, the same time series of sensor measurements would be obtained for two, three, and four injection locations since the latter two do not influence contamination measurements. If four sources are given to the prediction model as input, where contamination can be measured only from two sources, that would significantly reduce the accuracy of the prediction model. Thus, only nodes which contribute to the contamination measurements in sensors are considered for multiple injection scenario.

**Figure 1.** Net3 network with sensor layouts.

**Figure 2.** Richmond network detail with sensor layout.

An example of the proposed methodology can be seen for arbitrarily chosen Net3 contamination scenario in Figure 3. Randomly chosen contamination scenario parameters are 3 source nodes (159, 151 and 123), with contamination value of 200 mg/L, starting time 13 h and 20 min and injection duration 2 h. Sensor measurements for chosen contamination scenario can be seen in Figure 4. It can be observed that for source node 151 contamination scenario remains undetected in all sensors placed in the water distribution network, thus for multiple sources scenario only source nodes 123 and 159 are further considered.

**Figure 3.** Contours of chemical for randomly chosen Net3 contamination scenario 90 min after injection starting time. Contamination from source node 151 remains undetected, so the source node is not included for multiple injections scenario.

**Figure 4.** Sensor measurements for Net3 contamination scenario with (**a**) injection node 159, (**b**) injection node 123, (**c**) injection nodes 123 and 159 and (**d**) contamination measurements in the sensor in node 181.

#### *2.2. Demand Uncertainty and Sensor Type*

To investigate demand uncertainty, for both Net3 and Richmond networks, for every network node first it was randomly chosen if demand will be altered or not. If node base demand was to be altered, the percentage from 0–5% is randomly chosen for each network node, to reduce or increase base demand by the chosen percentage, resulting in a random demand span of 10%. To further investigate influence of demand uncertainty, the percentage from 0–10% is randomly chosen to reduce or increase base demand, resulting in a random demand span of 20%. All network demand patterns were kept the same, only base demand was changed. This method was conducted for every contamination scenario, thus resulting in different hydraulic conditions for each contamination scenario.

For sensor type influence, fuzzy sensor measurements were made where sensor detection was considered either low, medium, or high. Chemical concentration value *C* in range 0 < *C* < 300 mg/L was considered low, in range 300 < *C* < 1000 mg/L was considered medium and high if *C* > 1000 mg/L. Prediction model input features were defined as 0 if no contaminant was detected, 1 for low measurements, 2 and 3 for medium and high measurements, respectively.

## *2.3. Machine Learning Classifiers*

Two different machine learning classifiers, Random Forest and Artificial Neural Network were used to compare the efficiency of the proposed method. Random Forest al-

gorithm [28], based on multiple decision trees is used, with 250 estimators (trees) with a maximum depth of 30 and the minimum number of samples required to split an internal node 8. An artificial neural network with three hidden layers with 100 nodes in each layer, with hyperbolic tangent activation function and Adam solver for weight optimization is used. Proposed parameters were chosen with the grid search hyperparameter optimization method, while other parameters, which are not mentioned, are kept constant. Implementation in the Python library Scikit-learn [29] version 0.20.3 is used for both classifiers. Obtained data was split 70% for teaching and 30% for model testing. Flowchart of the prediction model can be seen in Figure 5. Data generation and prediction model training was done using the supercomputing resources at the Center for Advanced Computing and Modelling, University of Rijeka.

**Figure 5.** Flowchart of Machine Learning algorithm for prediction of number of contamination sources.

Input data for the prediction model is the time series of sensor measurements. For both Net3 and Richmond network, 25 features per sensor are obtained, which resulted in 100 features for Net3 and 125 features for Richmond network. The output of the machine learning model is the number of injection locations where two different prediction models are used. The first prediction model was used to predict the exact number of injection locations, i.e., 4 different classes are predicted. In the second model it is predicted only if single or multiple injections occurred, i.e., 2, 3 and 4 injection locations are treated as same, multiple injections class, thus only 2 different classes are predicted (single and multiple injections). To further increase the accuracy of the latter prediction model, the threshold value is introduced. Only if the model predicts a single source scenario with a probability greater than the chosen threshold value, single source prediction is made. In other cases, the scenario is treated as multiple sources. Threshold values of 50%, 60%, 70%, 80%, 90%, and 95% are investigated.

## **3. Results**

#### *3.1. Model Accuracy*

The influence of input data on prediction model accuracy is investigated for both benchmark networks where data ranged from 50,000 to 500,000 inputs (Figure 6). An investigation is conducted for prediction model with 2 categories (model predicts only if single or multiple injection locations are present) and with 4 categories (model predicts an exact number of injection locations). For each model and each number of inputs, 20 runs were conducted to take into consideration the influence of random seed. For the Net3 network second sensor layout with sensors placed in nodes 115, 119, 187, and 209 was considered. For Net 3 results are presented for both RF and NN prediction models. Standard deviation ranged from 0.63% for 50,000 to 0.33% for 500,000 inputs for NN model, and from 0.33% for 50,000 to 0.1% for 500,000 inputs. It can be observed that the RF model has slightly better accuracy for all investigated models. Also, due to the faster execution time of the RF model, for all further analyses, only RF results will be presented. For Richmond network, standard deviation ranged from 0.28% for 50,000 inputs to 0.12% for 500,000 inputs which indicates the stability of the model. Presented results are an average of all 20 runs.

**Figure 6.** *Cont.*

**Figure 6.** Accuracy of prediction models for different number of inputs for (**a**) Net3 network and (**b**) Richmond network.

It can be observed that even for a small number of input data considerable accuracy can be achieved. For model with 2 categories even with 50,000 inputs accuracy of the model is above 85% for both considered networks. After 200,000 inputs accuracy of the models for both networks tend to only slightly increase with the further increase of the number of input data. For 500,000 inputs accuracy of the Net3 network is 66.83% and for Richmond network 72.96%. When simplification is made, and the model only needs to predict single or multiple injection locations, accuracy significantly increases and for 500,000 inputs for the Net3 network is 91.46% and for the Richmond network 93.4%.

## *3.2. Threshold Influence*

To further increase the accuracy of the prediction model, the threshold value is introduced for the model which predicts 2 categories. Detailed results are presented for models with 500,000 inputs for Net3 (Tables 1 and 2) and Richmond network (Tables 3 and 4). Presented results are the average of values obtained from 20 runs. As expected, with the increase in threshold value accuracy of the prediction model increases. However, with a greater threshold value, a greater number of single injection scenarios, as a precaution, are classified as multiple sources, thus a smaller number of true single injection scenarios are detected. For both networks, when the threshold value is 95%, a very low percentage of correct prediction of single source scenarios can be observed when prediction model parameters chosen with grid search optimization method (250 estimators, maximum depth 30, minimum samples for split 8) were used (Tables 1 and 3). Thus, different prediction model parameters (180 estimators, maximum depth 80, minimum samples for split 10) were also investigated to test its influence on model accuracy when threshold values are considered. In Tables 2 and 4 it can be observed that for the greatest threshold value (95%) correct prediction of single sources scenarios greatly increases, and is around 30% of the total number of single source scenarios. As threshold value decreases, similar percentages are observed for both models, which indicates that model accuracy is similar for different RF parameters. However, when greater prediction certainty is expected, model parameters must be carefully considered.

For both networks, accuracy with threshold value 95% is above 99.5%. It can be observed from Table 2 that for Net3 only 36% of total number of single source scenarios are correctly predicted where for Richmond network (Table 4) that value is 37%. For threshold value 50% for Net3 94.5% of single injection scenarios are correctly predicted; however, the number of wrong predictions increases. The same can be observed for the Richmond network where for threshold value 50%, 97.8% of single injection scenarios are correctly predicted but the percentage of wrong single injection scenarios increases from 0.8% to 12.7%.

The problem remains with scenarios that are wrongly predicted even for a threshold value of 95%. With further increase of threshold value, the number of wrongly predicted scenarios would decrease, but only because ultimately all scenarios would be classified as multiple sources (this can also be observed in Tables 1 and 3 for first chosen RF parameters). Thus, optimum threshold value should be chosen to both provide a reasonable number of single injection scenario predictions but with a high model accuracy. In-depth analysis of scenarios where the model wrongly predicts a single injection scenario with a high threshold value should be conducted. Also, it should be investigated how much accuracy of the model can be further increased with a larger number of inputs and with the usage of different classifiers.

**Table 1.** Influence of threshold value on model accuracy for Net3 network (250 estimators, maximum depth 30, minimum samples for split 8). Percentage indicates number of predicted simulations based on total number of single source scenarios.


**Table 2.** Influence of threshold value on model accuracy for Net3 network (180 estimators, maximum depth 80, minimum samples for split 10). Percentage indicates number of predicted simulations based on total number of single source scenarios.


**Table 3.** Influence of threshold value on model accuracy for Richmond network (250 estimators, maximum depth 30, minimum samples for split 8). Percentage indicates number of predicted simulations based on total number of single source scenarios.



**Table 4.** Influence of threshold value on model accuracy for Richmond network (180 estimators, maximum depth 80, minimum samples for split 10). Percentage indicates number of predicted simulations based on total number of single source scenarios.

## *3.3. Sensor Layout*

The influence of sensor layout was tested for both Net3 and Richmond networks. 20 runs were conducted for the model with 500,000 inputs and average accuracy for all runs can be seen in Table 5. It can be observed that for the same number of sensors, their layout influences the accuracy of prediction models. This is expected, since the same behavior can be seen when the detection rate of contamination event is investigated for different sensor layouts. In the paper by Ostfeld et al. [30] for the same network and the same number of sensors detection likelihood of contamination event greatly differs for different sensor layouts. Results show that the prediction model for 2 categories (predicts single or multiple injections) is less influenced by sensor layout and all sensor layouts have accuracy around 90% or higher.

Interestingly, greater model accuracy can be observed when a smaller number of sensors is placed for Net3 layout with sensors in nodes 117, 143, 181, and 213 and for Richmond network. However, it can be explained with the fact that a greater number of contamination events remain undetected. i.e., with the greater number of sensors, contamination events from the greater number of network nodes are detected, resulting in more combinations when considering multiple injection locations. When sensor placement is sparser, a smaller number of network nodes can be detected when the contamination event occurs, resulting in a smaller number of combinations for multiple injection locations and consequently providing better model accuracy with 500,000 inputs.


**Table 5.** Influence of sensor layout for Net3 and Richmond networks on prediction model accuracy.

## *3.4. Demand Uncertainty and Fuzzy Sensors*

Influence of demand uncertainty and fuzzy sensors was investigated for Net3 network with 4 sensors in nodes 117, 143, 181 and 213 and for Richmond network with 5 sensors in nodes 93, 352, 428, 600 and 672. 20 runs were conducted for RF models with 500,000 inputs and average accuracy can be observed in Table 6. When demand uncertainty is considered the accuracy of RF models slightly decreases for both networks. The influence of fuzzy sensors is more prominent, where the greater reduction in prediction accuracy can be observed for the Net3 network. When considering both demand uncertainty and fuzzy sensors in the same model, accuracy further slightly decreases. However, it can be observed that for both networks model which predicts 2 categories has accuracy above 90% for all cases. This shows that the proposed model could be applied in a real case scenario.

**Table 6.** Influence of demand uncertainty and fuzzy sensors for Net3 and Richmond network on prediction model accuracy.


#### **4. Discussion**

Accuracy of prediction models for both networks has similar results with small differences, which shows that the proposed methodology could be successfully applied to other networks. Further investigation should be conducted for large size water distribution networks and different sensor placements, to fully investigate the robustness of the proposed method. Also, it must be noted that simplification was used in this study, where all source nodes had the same parameters (injection starting time, duration, and concentration value), thus, it should be investigated how the model predicts if those parameters are different for each injection node.

Although slightly, with the increase of input data model accuracy still increases, so in further study a greater number of data inputs should be investigated. Also, in the proposed scenarios report time step was chosen to be 1 h, resulting in 25 features per sensor. It should be investigated if a greater number of features, i.e., smaller report time step would increase model accuracy and if similar model accuracy could be achieved with a smaller number of contamination readings. The optimal number of features and inputs should be investigated to achieve great accuracy but with reasonable execution time. However, to obtain a greater number of inputs a greater amount of time is needed, so the model should be trained before the actual contamination event occurs. In that case, the model would be trained with simulation results with average demand patterns. This surely would mean that true contamination event will have different demands which would influence the accuracy of the prediction model. Investigation of demand uncertainty with arbitrarily chosen demand variation spans showed that small differences of base demands slightly influence prediction model accuracy. However, it must be taken into consideration that when base demand variation is defined with percentage, small demand variation is achieved when base demand is small and greater demand variation only when base demand is greater. Greater difference in demands should be further investigated since the usual variability of consumption can be greater than considered in this paper. Different machine learning models, with different expected demand patterns, can be prepared for contamination event so prediction can be obtained instantaneously. However, in case of contamination event, greater oscillations in the hydraulics of water distribution network could occur, such

as pipe burst or some other unplanned event, which would greatly influence change in demand patterns. Thus, it would be beneficial to investigate other algorithms that could increase accuracy with a smaller number of input data. In that case, input data can be obtained after the contamination event occurred, in a reasonable amount of time. That would be greatly beneficial since the simulation model can then be calibrated with sensor measurements from the field and input data would be more precise. The proposed method can be easily coupled with other machine learning approaches since inputs obtained for this model can also be used for teaching model that predicts injection location.

Investigation of different sensor layouts, demand uncertainty, and fuzzy sensors showed that sensor layout and type of sensors have the greatest impact on prediction model accuracy. Demand uncertainty slightly decreases model accuracy. However, model accuracy can be greatly reduced when a real case event is considered since both demand uncertainty and measurement errors can be greater than considered in this work. Thus, a threshold value is introduced which can help increase model accuracy. Greater threshold value increases model accuracy; however, it also leads to a greater number of single injection scenarios classified as multiple injections. It is also observed that prediction models are not very sensitive to model parameters; however, when threshold value is used, i.e., model prediction certainty is evaluated, model parameters are very important for method efficiency. Thus, the investigation of different machine learning approaches should be further investigated to increase model accuracy.

When observing presented results it must be taken into consideration that numerical model simplifications are made, where EPANET was used which assumes complete mixing in all network junctions and uses pure advection transport model. Also, in the presented study benchmark networks are used, and numerical simulations are conducted for only 24 h, where more than 24 h are needed to obtain stable contamination scenario results. However, the functionality of the presented machine learning approach is not dependent on the numerical model setup, and it is assumed that the same numerical approach that is chosen for the optimization process is to be also chosen for the prediction model preparation. In this way, all discrepancies due to numerical model simplifications would be also present in the optimization and as such are not the result of using the proposed machine learning approach. Furthermore, network uncertainties were not considered regarding internal pipe diameter and pipe roughness which should be considered in the further research.

## **5. Conclusions**

In this paper, the machine learning approach is presented which helps identify the number of injection locations based on sensor measurements. Random Forest classifier and Neural Network classifier are used on medium-sized benchmark network, where Random Forest classifier provided better accuracy and faster execution time, thus is used for all other investigations. Two different sized benchmark networks are considered, where it is shown that the machine learning approach can be successfully used to predict the number of injection locations. This can help define the number of optimization parameters, where redundant parameters can be avoided which needlessly increase the complexity of the problem. The prediction model shows great accuracy when it predicts only if single or multiple injection locations occurred. The threshold value is proposed which further increases model accuracy since the single injection scenario is assumed only if the model predicts with certainty greater than the threshold value. Lower accuracy is obtained when the exact number of injection locations is predicted. The accuracy of the prediction model is investigated for different sensor layouts and in case of demand uncertainties and fuzzy sensors. Conducted research showed promising results, where exploration of other algorithms and increased number of input data should be investigated to further increase the accuracy of both models.

**Author Contributions:** Conceptualization, I.L. and L.G.; Data curation, I.L.; Formal analysis, I.L; Investigation, I.L. and L.G.; Methodology, I.L. and L.G.; Resources, Z.C. and L.K.; Software, I.L.; ˇ Supervision Z.C. and L.K.; Validation, I.L.; Visualization; I.L.; Writing—original draft, I.L.; Writing— ˇ review and editing, L.G., Z.C. and L.K. All authors have read and agreed to the published version of ˇ the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **A Machine Learning-Based Algorithm for Water Network Contamination Source Localization**

## **Luka Grbˇci´c 1,2, Ivana Luˇcin 1,2, Lado Kranjˇcevi´c 1,2,\* and Siniša Družeta 1,2**


Received: 3 April 2020; Accepted: 30 April 2020; Published: 3 May 2020

**Abstract:** In this paper, a novel machine learning based algorithm for water supply pollution source identification is presented built specifically for high performance parallel systems. The algorithm utilizes the combination of Artificial Neural Networks for classification of the pollution source with Random Forests for regression analysis to determine significant variables of a contamination event such as start time, end time and contaminant chemical concentration. The algorithm is based on performing Monte Carlo water quality and hydraulic simulations in parallel, recording data with sensors placed within a water supply network and selecting a most probable pollution source based on a tournament style selection between suspect nodes in a network with mentioned machine learning methods. The novel algorithmic framework is tested on a small (92 nodes) and medium sized (865 nodes) water supply sensor network benchmarks with a set contamination event start time, end time and chemical concentration. Out of the 30 runs, the true source node was the finalist of the algorithm's tournament style selection for 30/30 runs for the small network, and 29/30 runs for the medium sized network. For all the 30 runs on the small sensor network, the true contamination event scenario start time, end time and chemical concentration was set as 14:20, 20:20 and 813.7 mg/L, respectively. The root mean square errors for all 30 algorithm runs for the three variables were 48 min, 4.38 min and 18.06 mg/L. For the 29 successful medium sized network runs the start time was 06:50, end time 07:40 and chemical concentration of 837 mg/L and the root mean square errors were 6.06 min, 12.36 min and 299.84 mg/L. The algorithmic framework successfully narrows down the potential sources of contamination leading to a pollution source identification, start and ending time of the event and the contaminant chemical concentration.

**Keywords:** machine learning; artificial neural networks; random forests; water network pollution; sensor networks; parallel computing

## **1. Introduction**

Identifying the source of contamination in a water supply network is an important task since a contamination event is potentially hazardous to the human population in an urban environment. Additionally, a fast identification of a pollution source enables the governing authorities to rapidly react in order to stop the further spread of the contaminant through the water supply network.

Researches have tackled the issue of water supply networks contamination in several ways which include an optimal positioning of water quality sensors in a network [1–3] to facilitate the source identification process and optimally cover all possible intrusion points, rapid contamination event response procedures [4–6] and simulation-optimization methods for contamination source detection and duration based on simulation of the water network contamination event [7–9]. Many researches have incorporated additional uncertainties into the hydraulic simulation process which include uncertain sensor measurements and water demand variability [10–12]. A recent and thorough

review of various approaches in tackling the water supply pollution source identification problem can be found in Adedoja et al. [13].

The optimization-simulation approach for finding the source of pollution in a network entails that an optimization algorithm is being coupled with a water supply network hydraulic simulator and the difference between the measured sensor data and the optimization algorithm generated and simulated values (sensor water quality readings through a certain time interval) is being minimized. Through this procedure, the source of contamination, starting time and the end time of the contamination event and the contaminant chemical concentration are obtained.

Probabilistic approaches are also possible for determining the source of contamination in a water supply network and the methods for this approach that were used in previous studies include the Bayesian Belief Networks (BBN) [14–16] and backward probabilistic models [17] which show the ability to predict the source of the contamination with a high probability. Data-driven methods are also a possible tool for the water supply network contamination problem. In [18,19], a database was compiled by massive data mining of hydraulic and water quality simulation contamination events for a fast identification of a source in case of a real event. The sensor readings in case of a contamination event are matched with the simulation events from a database and in [19] a statistical maximum likelihood approach was used for matching. In [20], Monte Carlo (MC) water quality and hydraulic simulations were run in parallel and then used to detect the source of pollution with a certain criteria.

In [21] the logistic regression approach was used to determine top candidates for being the true contamination source with them additionally being explored with local search methods to determine other relevant variables of a contamination event (start time, end time and chemical concentration). The input data for the logistic regression were the sensor readings through time that were constantly updated with new data.

Data-driven models would include using machine learning algorithms to localize the contamination sources in a water supply network. Kim et al. [22] used an Artificial Neural Network (ANN) to find the source of pollution in a small network and Rutkowski and Prokopiuk [23] used a learning vector quantization Neural Network (LVQNN) to locate a zone with a supply network where a potential source of contamination would be located. Wang et al. [24] used Least Squares Support Vector Machines (LS-SVM) to enhance the reliability and accuracy of water sensor contamination detection.

Previous studies have extensively explored Monte-Carlo based methods in air pollution and groundwater pollution source detection problems. In Guo et al. [25] a Markov Chain Monte Carlo (MCMC) sampling method coupled with a Bayesian probabilistic approach was applied to find a source of unsteady atmospheric dispersion which was numerically modeled. Wade and Senocak [26] used Bayesian inference MCMC for the purpose of reconstructing multiple air pollution sources. The study was tested on real-field data and the method includes a ranking of the most probable number of pollution sources which is based on error analysis. In the work by Bashi-Azghadi et al. [27] Probabilistic Support Vector Machines (PSVM) and Probabilistic Neural Networks (PNN) were used to determine an unknown source of pollution in a groundwater system. Vesselinov et al. [28] studied the application of semi-supervised machine learning methods with synthetic and real measured data to identify contamination sources of chemical mixtures in groundwater flows.

In this study we present an algorithm which utilizes the ANN for classification of contamination source in a water supply network and Random Forests (RF) for prediction of contamination start time, end time and chemical concentration. The algorithm is built in a high performance computing (HPC) environment and uses a parallel tournament style selection of most probable contamination source node between a group of nodes. All network nodes are divided in chosen number of tournament groups, where each group is assigned to a single processing core. Monte Carlo (MC) hydraulic and water quality simulations using EPANET2 (Rossman [29]) are run in parallel and the obtained simulation results are used to create models for every tournament group. Network nodes are randomly distributed into tournament groups and the parameters for every simulation (contamination start time, end time and chemical concentration) are randomized. Each network node in a tournament

group obtains the same number of results (or MC simulations). Once a node is selected as a winner in the tournament process from existing groups, new tournament groups are created with a reduced number of suspect nodes and this process is repeated until a stopping criterion is satisfied which in this case is the number of set tournament loops which is also a parameter of the algorithm. The whole algorithmic framework was built with a combination of the Python 3.7 programming language and Simple Linux Utility for Resource Management (SLURM) for HPC systems. The ANN classification and the RF regression analysis were done using the Python machine learning library scikit-learn 0.22. The algorithm was tested on two benchmark networks taken from previous studies and it shows good results but also a possibility for improvement with future studies. The algorithm shows a great ability to include the true source node as a top candidate among the remaining source nodes and is good at predicting the other relevant contamination event parameters which can be further used for coupling with optimization algorithms.

## **2. Materials and Methods**

## *2.1. Water Supply Network Benchmarks*

The algorithmic framework was tested on two benchmark networks—the Net3 example from EPANET2 and the Richmond water supply network (Van Zyl [30]).

The Net3 EPANET2 benchmark water supply network consists of 92 nodes and was specifically made for water quality hydraulic simulations. The simulation parameters are set as: the total simulation time is 24 h with a 10 min hydraulic time step and 5 min water quality time step and a 10 min pattern time step. The already optimized water quality sensor layout was set as the one from the work by [7] (network nodes 117, 143, 181 and 213 are set as sensors). The Net3 network layout with the sensor placement can be seen in Figure 1. In each MC simulation of Net3, both the start and the end times (*Sm*, *Em*) of the contamination event were randomly set from 0 to 24 h with an obvious restriction that *Em* > *Sm*. The value of *Cm* was randomly chosen from an interval from 10 to 1000 mg/L and was kept constant throughout the whole contamination scenario. The sensors used in Net3 recorded data during the whole 0–24 h interval for every hour (a total of 25 water quality measurements per sensor) which means that there were 100 input features for the RF regression analysis.

The Richmond network consists of 865 nodes and it was obtained from The Centre for Water Systems (CWS) at the University of Exeter [31]. Simulation time was set as 72 h with a 1 h hydraulic time step, a 5 water quality time step and a 1 pattern time step. The sensor layout was set according to the work by Preis and Ostfeld [7] (nodes 123, 219, 305, 393 and 589 are sensors nodes). The Richmond water supply network layout detail can be seen in Figure 2. The selection of random parameters for each MC Richmond network simulation were defined the same way as for the Net3 network. The Richmond network sensors recorded data during the 0–72 h interval for every hour (a total of 73 measurements per sensor) making a total of 365 input features for the RF regression analysis.

**Figure 1.** Net3 water supply network layout with sensor positioning by Preis and Ostfeld [7].

**Figure 2.** Richmond water supply network layout with sensor positioning by Preis and Ostfeld [7].

## *2.2. Algorithmic Framework*

The algorithm presented in this study was designed to work in a HPC environment in order to detect the water supply network contamination source in a rapid and efficient way. It is based on distributing all potential contamination source nodes of a water supply network into subgroups for which MC simulations would be done in conjunction with machine learning methods to determine which node would best fit to be the true source node of contamination. In this way, each node in a subgroup would be a part of a tournament in which the node with the highest probability of being the contamination source would continue to the next tournament round. After each tournament, the winning nodes would be redistributed in a new tournament group until a predetermined number of tournaments was reached. Each CPU in a HPC system is assigned a tournament node group and after each tournament the total number of used CPUs would decrease since the losing nodes would be discarded and the remaining nodes which are fewer would form new tournament groups. The flowchart of the whole algorithmic framework is shown in Figure 3.

**Figure 3.** Algorithm flowchart.

As seen in the flowchart, the algorithm is initialized with reading all potential source nodes X in a water supply network and then distributing them in tournament groups of constant size k which is a parameter of the algorithm and can be freely selected. As each tournament group is assigned to a CPU, the number of used CPUs n is determined with n = X/k. After distributing the tournament groups to each CPU, MC simulations are performed m times (m/k times for each suspect node in the tournament group with ideally the modulus of m/k being zero—if not, a node is randomly selected for the additional run) with randomly selected starting *Sm* and ending *Em* times of the contamination event, and the contaminant chemical concentration *Cm*. Input and output data of each MC simulation are being saved for each suspect node x in every tournament group n.

After the MC simulations are done, the input (sensor readings through time) and output data (source node with used *Sm*, *Em* and *Cm*) are used for training each tournament group's machine learning (ML) model. The ML model used can be any ML classifying algorithm which supports the prediction of multiple classes. ML output variable set consists of all network nodes that are within that tournament group. After the model was trained with the MC generated data, the sensor readings of the contamination event are being used for the ML prediction of the most probable source node in each tournament group. The nodes with the highest probability in all tournament groups are considered to be the tournament winners. After every used CPU generated a tournament winner, a list of all winners is compiled and if the number of set tournament loops l is not exceeded, the tournament process and distribution is repeated and the number of nodes X is updated (it is equal to the number of winners). In this case, the number of X should be smaller than in the previous tournament loop and consequently the number of used CPUs is reduced since it is dependant on the number of nodes for distribution. It is important to note that each winning node's input and output data is saved from every tournament loop.

If the freely selected algorithm parameter l is exceeded, each winning node is then assigned again to a CPU and with its previously obtained MC input and output data, a ML model is trained and a prediction is performed for the remaining variables of the contamination event (*Sm*, *Em* and *Cm*). The predicted values of *Sm*, *Em* and *Cm* of each winning node's ML regression model are then used for simulating the contamination event scenario and the obtained sensor readings are then compared with the real contamination event sensor readings with a RMSE analysis, which in turn creates a ranking where the node with the smallest RMSE is placed at the top.

This whole algorithmic framework was built within the open source workload manager for cluster systems SLURM and the Python 3.7 programming language.

## *2.3. ANN Classifier*

In the previous sub section, the tournament ML classifier was generally defined in the whole algorithmic framework and basically any ML algorithm which can predict multiple classes can be used. In this study, the ANN algorithm is used for classifying the most probable source nodes in a tournament group.

The Multi-layer Perceptron (MLP) type of ANN was used from the Python 3.7 machine learning library scikit-learn 0.22 [32]. The MLP ANN was constructed with both input and output layers and three hidden layers in between. Both first and last hidden layers consisted of 100 neurons, while the middle hidden layer was formed with 500 neurons. The stochastic gradient-based optimizer ADAM for MLP weights optimization was selected through the process of hyper parameter tuning, just as the number of neurons in every ANN layer.

With a preliminary analysis of the possible input variables of the ANN MLP model it was determined that great accuracy of the model can be achieved if only the maximum values of the chemical concentration recorded per sensors through a time interval in the water supply network are used. The preliminary analysis was done through 10 runs and each run when the true source node was in the top 6 of the final nodes was considered successful. The analysis was done on the Net3 benchmark network with the contaminant source node 119 as described in Section 3.1. The true source node was a part of the top 6 ranking suspect nodes for all of the 10 preliminary runs.

This means that the number of neurons at the input layer is equal to the number of sensors used in the water supply network. Furthermore, using the whole time interval of all sensor water quality readings as ANN MLP inputs was tested in the preliminary analysis and it was found that the performance was not better (8 out of 10 runs were successful) than using the maximum values of recorded water quality through time, so naturally, the maximum values per sensor were used as inputs since the number of ML model features is much smaller that way.

The output of the MLP was a list of all tournament group nodes and an assigned probability of each node being the true contamination source after the real contamination event sensor reading was evaluated with the trained model.

In Figure 4 the whole MLP can be seen with *max*(*Csn*(*t*)) being the maximum concentration recorded by the n-th sensor in the network and *Nn*% being the probability that the tournament node n is the true source node. All of the MC generated input and output data (of each tournament group) is used to train the ANN model for each tournament group as the goal of the classifier is to predict the most probable contamination source node of each tournament group. The success was assessed by observing the prediction of the whole algorithmic process and not the accuracy of each tournament group ANN model.

**Figure 4.** MLP with all layers.

## *2.4. RF Regression*

The ML regression model for each tournament winning node after the parameter l was exceeded in the algorithmic framework can also be done with any ML algorithm which supports multi output regression. The RF algorithm (Breiman [33]) from scikit-learn 0.22 was selected for this purpose. All parameter values of the algorithm were set as default except for the number of estimators (trees in the random forest) which was set to be 200 with the process of hyperparameter optimization.

The input values for the RF regression analysis were sensor water quality readings throughout the whole time interval of the simulation (unlike the inputs used for the MLP ANN) and the output variables were the predicted values of *Sm*, *Em* and *Cm* for every winner node. A flowchart of the RF regression is seen in Figure 5 with *Csn*(*t*0...*tmax*) representing the chemical concentration recorded by the sensor n during simulation time.

**Figure 5.** Flowchart of the RF regression analysis.

## **3. Results and Discussion**

## *3.1. Net3 Network Contamination Scenario*

The contamination event scenario for the Net3 benchmark network was chosen to be from the same node (119) as in the one from the work by Preis and Ostfeld [7] and the location can be seen in Figure 6. The contamination event characteristics at source node 119 were freely chosen with the event starting at 14:20 h and lasting until 20:20 h with a constant chemical mass inflow of 813.7 mg/L.

The selected number of algorithm loops l was 3, the number of m (MC simulations for every tournament group) was 200 and the size of a tournament group k was 2, which means that with 92 initial water supply network nodes, the number of used CPUs for every tournament group was 46 and after every loop that number was halved. After three loops, the number of tournament winners was 11, which means that 11 CPUs were used for the RF regression analysis and prediction of other relevant variables.

**Figure 6.** Net3 water supply network contamination source location.

The contamination source search was repeated 30 times since there the algorithm consists of a stochastic component (MC simulations). In 22 out of 30 runs the true source node was the suspect node with the highest probability and in the remaining eight runs the true source node was a part of the final winners list, which means that the ANN classification can successfully narrow down the search space from 92 to 11 nodes in this case. The average run which includes MC simulations, ANN classification and RF regression lasted for 8 min (even though the RF regression lasted only for 8 s). The algorithm was run (at its initial loop) on 46 Intel Xeon E5 CPUs (two cluster nodes). Out of the other eight runs when the true source node was not ranked first, it was always in the top six of the tournament winners.

In Figure 7 a comparison can be seen between the true contamination event (14:20 h to 20:20 h with 813.7 mg/L) and all of the 30 predicted contamination events for the true source node. It can be seen that the end time prediction of the event is very accurate while the starting time only lacking in accuracy on three runs. The overall RMSE for the starting time for all 30 runs is 48 min, the end time is 4.38 min and the chemical concentration is 18.06 mg/L. The average RMSE For the three of the worst runs with respect to the starting time was 2.47 h. In Table 1, a summary of all runs can be seen through the RMSE analysis and the successful runs represent how many times of the total of 30 runs the true source node was part of the final tournament. The minimum and maximum errors for *Sm*, *Em* and *Cm* for all 30 runs are presented in Table 2.

**Figure 7.** Net3 network true contamination event (black) and predicted contamination events (grey).

In Table 3 the best and worst runs are compared with the true contamination event parameters for *Sm*, *Em* and *Cm*. The overall best and worst runs are calculated (individual RMSE) by taking into account all of the three variables.

**Table 1.** Contamination event RMSE analysis for Net3 network of all 30 runs.


**Table 2.** Net3 network minimum and maximum errors (*Sm*, *Em* and *Cm*) for 30 runs.


**Table 3.** Net3 network contamination event results comparison between true, best and worst of the total 30 runs.


In Figure 8 the nodes which were ranked first in the 30 runs can be seen along with a corresponding number of times they were ranked first. It can be observed that the nodes are topologically clustered together. This is expected since due to the multimodal nature of the problem.

**Figure 8.** Number of times out of 30 runs for which a node (marked blue) was ranked as first for the Net3 contamination event.

## *3.2. Richmond Network Contamination Scenario*

The Richmond network contamination event scenario was chosen to start at the same node (153) as in the work by Preis and Ostfeld [7] and the location can be seen in Figure 9. The contamination event characteristics at source node 153 were chosen with the event starting at 06:50 h and lasting until 07:40 h with a constant chemical mass inflow of 837 mg/L.

The selected number of algorithm loops l was 1, the number of m (MC simulations for every tournament group) was 2500 and the size of a tournament group k was 4, which means that with 865 initial water supply network nodes, the number of used CPUs for every tournament group was 432. After 1 loop the number of tournament winners was 217, which means that 217 CPUs were used for the RF regression analysis and prediction of other relevant variables.

**Figure 9.** Richmond water supply network contamination source location.

The contamination source search was repeated 30 times just like for the Net3 contamination search. The true source node was ranked first in seven out of 30 runs and it was in the tournament winners list 29 out of 30 times which means that the true source node was not subjected to the RF regression analysis for only one run. The total average run time of the MC simulations and the ANN classification was 41 min and 75 s for the RF regression analysis. The algorithm was run (at its initial loop) on 432 Intel Xeon E5 CPUs. Even though it was ranked first in only seven runs, that was the most number of times a node was ranked first out of the 30 runs. In the 29/30 runs it was in the winners list (which consisted of 217 nodes) it always finished in the top 10 after RF regression was completed.

In Figure 10 the comparison between the true contamination event and the predicted contamination events can be seen. The starting and ending times RMSE for the 29 of the 30 total runs are 6.06 min and 12.36 min respectively, while the chemical concentration RMSE is 299.84 mg/L. The starting and ending times of the predictions are in good agreement with the true values but the chemical concentration value is underestimated by the RF regression.

A RMSE analysis summary of all runs with the number of successful runs can be seen in Table 4. The minimum and maximum errors for all 29 runs for the Richmond network are shown in Table 5 and in Table 6, the best and worst runs comparison is shown in terms of all the predicted variables (*Sm*, *Em* and *Cm*).


**Table 4.** Contamination event RMSE analysis for Richmond network of all 30 runs.

**Table 5.** Richmond network minimum and maximum errors (*Sm*, *Em* and *Cm*) for 30 runs.


**Table 6.** Richmond network contamination event results comparison between true, best and worst of the total 30 runs.


**Figure 10.** Richmond network true contamination event (black) and predicted contamination events (grey).

## *3.3. Algorithm Parameters Investigation*

In this subsection, an analysis of the influence of algorithm parameters is given. The required number of MC simulations m for each node, the tournament group size k and the number of algorithm loops l are separately explored. The examination of m, k and l was done on the Net3 water supply network with the previously defined contamination event scenario (source at node 119, *Sm* = 14:20 h, *Em* = 20:20 h and *Cm* = 813.7 mg/L). Each different setup (of m, k and l) was independently run 10 times and a run was considered successful if the node was present in the final ranking after the RF regression and RMSE analysis.

Firstly, the influence of parameter m on the prediction accuracy was explored and a complete summary can be seen in Table 7. The selected tournament group size k for all runs was 2 and the number of loops l was set as 1 during the exploration of parameter m.


**Table 7.** Influence of the number of MC simulations on the accuracy and efficiency of the algorithm.

The first column of Table 7 represents the total number of MC simulations per tournament group, which means that since the tournament group size was 2 for each run, each node in the tournament

group was the source node in m/2 MC simulations. When observing the analysis of the parameter m in Table 7 it can be seen that the higher the value of the parameter is, both contamination event prediction accuracy (as seen through the decrease *Sm*, *Em* and *Cm* RMSE) and the average computation time per run (last column) increases. This is expected since more randomly generated data covers more possible scenarios and more input data for the ML model enables a wider and more accurate solution space exploration. Additionally, the best and worst ranks of the true source node are shown and the number of times the true source node won (Times won meaning the rank was 1).

When m was set as 80 the number of successful runs was 10 and the worst possible rank was 6. This means that after 2 min of computation time on average, the search space was reduced from 92 total nodes to 6, which is a reduction of 93.5%. The set of runs when m was 400 could be considered as the first set of runs when the results are acceptable in terms of finding the true source node since it was the winner 50% of the time.

When the value of m is 2000 and above it can be seen that there is not a significant change in prediction accuracy as the RMSE values of *Sm* and *Cm* exhibit stability and minor oscillations, while *Em* has showed a steady convergence to the same value as the true contamination scenario for all 10 runs.

For further exploration of the tournament group size k, the chosen m was 800 as it exhibited a reasonable computation run time, accuracy in terms of the average RMSE and the number of times the true source node was ranked first. The same scenario was chosen as the one for the investigation of m with the number of tournament loops l set as 1. Five different tournament group sizes k were explored and are summarized in Table 8.

**Table 8.** Influence of the number of the tournament group size k on the accuracy and efficiency of the algorithm.


From the results presented in Table 8 it can be observed that when the tournament group is larger, both accuracy and prediction reliability decrease. Furthermore, besides tournament group sizes of 2 and 4, a reasonable result in terms of reliability is achieved with k = 10 with a total search space reduction of 94.6% and even a 50% winning rate in the 10 successful runs. The number of used CPUs for each tournament group size is added and with the given as the last column of the Table 8. Even though a tournament group size of 2 is not that impressive when compared to those of 4 and 10 in the categories of best and worst rank and times won, the achieved overall *Sm*, *Em* and *Cm* RMSE shows that it is undoubtedly more accurate.

Lastly, the influence of the number of tournament loops l is investigated and a summary of the results is shown in Table 9 . The same scenario was used as for the exploration of previous two parameters with k = 2 and with a total number of MC simulations m = 800. An additional column m/L was added to the Table 9 which defines the number of MC simulations m (of a tournament group) per every loop l.

It can be observed that increasing the number of loops l up to a certain value increases the accuracy and reliability of the algorithm. Even though the total number of MC simulations is the same for every run and the computational strain in that sense is similar, adding more loops decreases the number of used CPUs after every tournament loop since losing nodes are omitted and that can be considered as a great advantage.


**Table 9.** Influence of the number of the tournament loops l on the accuracy and efficiency of the algorithm.

The value of *Sm*, *Em* and *Cm* RMSE does not differ much for all tested loops l since the total number of MC simulations is preserved. When the number of loops is set to 8, the successful number of runs dropped as the number of MC simulations per loop (m/L) was not high enough and the true source node was not in the final RF analysis ranking for one run. This was also observed for smaller values of m in Table 7. It can be argued that a higher number of loops positively affects the success of the algorithm in predicting the relevant variables (source node, *Sm*, *Em* and *Cm*) since there is a higher chance that a main tournament top ranking rival to the true source node is omitted in the process of removing losing nodes after every tournament loop. However, setting l too high could result in unsuccessful runs as well (due to a small m/L) as it can be seen in Table 9.

#### **4. Conclusions**

In this paper a novel algorithmic framework for water supply network contamination source node identification is presented and tested on two different benchmark networks. The algorithm is specifically created for massively parallel HPC systems and it utilizes a combination of MC simulations and ML methods to identify the contamination event source node and all the relevant variables such as starting and ending times of the event and the contaminant chemical concentration.

The algorithm is based on running an equal number of MC simulations on a group of nodes in parallel and then selecting the node (via MLP ANN) with the highest probability of being the true contamination source node as the winner of the group (tournament). The number of MC simulations per tournament group is set as a parameter of the algorithm just like the number of tournament loops and the size of a tournament group. After a set of tournament winners is created, the algorithm utilizes a ML regression analysis using RF in parallel and creates a ranking of potential source nodes based on a simple error analysis with the true contamination event sensor data recorded.

The novel algorithmic framework, tested on two realistic and complex benchmark networks cases, displays the capability to narrow down the search space for the source node efficiently, leading to a pollution source identification. The algorithm can be also used in predicting the starting and ending times of the contamination event and the contaminant chemical concentration.

An investigation of algorithm's parameters which are the number of MC simulations, size of a tournament group and number of loops was also conducted. It is demonstrated that increasing the number of MC simulations is beneficial to the algorithm's ability to predict the true source node and relevant variables since more randomly generated data entails a broader solution space coverage, however this comes with an increase in computational time. It was observed that in order to increase the reliability of an accurate prediction the size of a tournament group should be as low as possible, depending on computational resources. Increasing the number of tournament loops shows to be advantageous in prediction accuracy; however, the number of MC simulations for each tournament loop should be high enough in order to preserve the reliability of prediction.

Further research is needed in determining the connection between the network size (initial number of potential source nodes) and the number of MC simulations needed to cover the search space of the contamination event efficiently and thoroughly. Furthermore, additional research should be done regarding the used ML methods investigating the use of different classifiers for the tournament group winner selection and various ML algorithms capable of multi-output regression analysis for the final node ranking. It would be also possible to couple the algorithm with simulation-optimization methods for an even faster convergence towards a true pollution source node detection in a way that only the tournament winners are subjected to the optimization process.

**Author Contributions:** Conceptualization, L.G., I.L. and S.D.; Data curation, L.G.; Formal analysis, L.G.; Funding acquisition, L.K.; Investigation, L.G. and I.L.; Methodology, L.G. and L.K.; Project administration, L.K.; Resources, L.K.; Software, L.G.; Supervision, S.D.; Validation, L.G.; Visualization, L.G.; Writing—original draft, L.G.; Writing—review & editing, I.L., L.K. and S.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Center for Advanced Computing and Modelling, University of Rijeka.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Abbreviations**

The following abbreviations are used in this manuscript:


## **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Extraction of Land Information, Future Landscape Changes and Seismic Hazard Assessment: A Case Study of Tabriz, Iran**

## **Ayub Mohammadi 1, Sadra Karimzadeh 1,2,3,\*, Khalil Valizadeh Kamran 1,\* and Masashi Matsuoka <sup>3</sup>**


Received: 7 October 2020; Accepted: 7 December 2020; Published: 8 December 2020

**Abstract:** Exact land cover inventory data should be extracted for future landscape prediction and seismic hazard assessment. This paper presents a comprehensive study towards the sustainable development of Tabriz City (NW Iran) including land cover change detection, future potential landscape, seismic hazard assessment and municipal performance evaluation. Landsat data using maximum likelihood (ML) and Markov chain algorithms were used to evaluate changes in land cover in the study area. The urbanization pattern taking place in the city was also studied via synthetic aperture radar (SAR) data of Sentinel-1 ground range detected (GRD) and single look complex (SLC). The age of buildings was extracted by using built-up areas of all classified maps. The logistic regression (LR) model was used for creating a seismic hazard assessment map. From the results, it can be concluded that the land cover (especially built-up areas) has seen considerable changes from 1989 to 2020. The overall accuracy (OA) values of the produced maps for the years 1989, 2005, 2011 and 2020 are 96%, 96%, 93% and 94%, respectively. The future potential landscape of the city showed that the land cover prediction by using the Markov chain model provided a promising finding. Four images of 1989, 2005, 2011 and 2020, were employed for built-up areas' land information trends, from which it was indicated that most of the built-up areas had been constructed before 2011. The seismic hazard assessment map indicated that municipal zones of 1 and 9 were the least susceptible areas to an earthquake; conversely, municipal zones of 4, 6, 7 and 8 were located in the most susceptible regions to an earthquake in the future. More findings showed that municipal zones 1 and 4 demonstrated the best and worst performance among all zones, respectively.

**Keywords:** remote sensing; GIS; Markov chain; land use; urban information; Tabriz City

## **1. Introduction**

Rapid urbanization, deforestation and increasing population have led to global environmental changes [1]. Because of this, large areas of agricultural land are being converted into urban land and industrial estates, which are prone to land degradation [1,2]. Typically, urbanization influences climate and water quality [3], which can result in changes in local climate. One of the main means by which to understand the relationship between humans and their environment is recording the changes occurring where they live [3,4]. Useful information can be obtained from the pattern and direction of land cover changes, from which better planning for sustainable development is possible. Landsat satellite imageries are widely used medium-scale data for land surface change analysis [2,5]. The low temporal baseline of these data is considered as a weak point that sometimes results in the omission of some dynamic changes [3]. Urban sprawl significantly changes landscapes across urban areas [6], which is usually associated with changes among vegetation, built-up areas and bare lands in the near or remote future. For seismic hazard assessment caused by an earthquake, urban information is urgently needed [7]. Normally, field checks and digital interpretation using the technology of remote sensing (RS) are among the common ways to extract urban information [1,7,8]. Geographic information systems (GIS) together with remote sensing technology can provide useful information for decision-makers [9]. RS data have successfully been applied for mapping and measuring the area and the extent of land cover. Satellite image performance has now been improved to a ground resolution of less than 1 m to be acquired. Physical assessment of urban areas from remotely sensed data enables comparative analysis of a city's extent within a region [3]. In order to record all changes, selecting the proper temporal baseline of RS data is very important [10]. The Landsat satellite imagery (medium spatial resolution RS data) has widely been utilized to quantitate urbanization across the world [2,11].

Predicted land use changes can help land use planners in mitigating the negative impacts on the environment. In recent years, the city and the surrounding areas have been several times shocked by large and small earthquakes. In developed countries, building inventory information is provided by local and central institutes [7]. Nevertheless, in developing countries like Iran, this is quite the opposite, so researchers should work and provide information to different organizations. Normally, this kind of study requires great effort and considerable financial support. There are many faults in Iran, a few of which have recently been activated and claimed many lives and also caused a great deal of damage to properties [12]; because of the possibility of the recurrence of such events in the near or remote future, fear still exists among inhabitants in these areas. Typically, because of the high buildings, this fear is most common among people in large cities like Tabriz.

Efforts have been made to carry out land cover information extraction using RS data and techniques. Urbanization influences ecosystems, but in order to determine and understand these impacts, precise and accurate information about the land cover's temporal and spatial changes is essential [2,3]. Built-up areas are very important for many studies, including those considering buildings' age, seismic hazard assessment, and prediction, so as to enable the optimal updating of built-up areas. Maximum likelihood, as one of the best methods for classification [11,13], was employed for extracting land cover using four cloud-free Landsat data. Accordingly, building inventory and urban sprawl information are important factors for damage estimation [7,14,15]. In order to determine if older buildings have been constructed under older seismic hazard standards, a map of the ages of buildings is needed [7,12]. For a fast and effective response during an earthquake, an urgent evaluation is needed [7], but, at present, sufficient information is provided to be used in the future to mitigate possible damages in the study area. Previous research simply studied the land cover detection of the city, but a comprehensive study on land cover change detection, future potential land cover, municipal performance evaluation and seismic hazard assessment is missing for Tabriz, which is one of the largest and most important cities in Iran. Therefore, this study presents a more comprehensive examination of the study area compared to the previous research, making it highly significant.

Many studies conducted land cover change detection and prediction using different models of fuzzy logic modeling [16–20]; geo-statistical methods [21–24]; Markov-CA [5,10,25–32]; cellular automata models [33–36]; propagating aleatory and epistemic uncertainty [37,38]; artificial neural network [39–43]; Hopfield neural network [44–48]; supervised back-propagation neural network [49,50]; self-adaptive cellular based deep learning [51–54]; analytical hierarchy process [27,55,56]; geographic information system (GIS)-based hybrid site condition [15,57,58]; recurrent neural network [59–62]; change vector analysis [63–67]; and different satellite imageries of both SAR and optical [1–3,6,7,9–11,13,14,68–76].

Previous studies on Tabriz City were not comprehensive, failing to include seismic hazard assessment and municipality performance. In recent years, the rapidly increasing population of the city has necessitated the construction of high-rise buildings and conversion of agricultural land into built-up areas. A comprehensive study that exploits the advantages of both remote sensing and GIS can turn satellite data into an actionable level that can be used for proper environmental planning. For the study area, a comprehensive evaluation from the point of view of satellite datasets and GIS is lacking. To fill this gap, GIS layers of various spatial information, SAR images from Sentinel-1 and optical images from Landsat missions were gathered to perform a new study of land cover change detection, future potential landscape, distribution of buildings by age, municipal performance evaluation and building damage assessment in Tabriz.

## **2. Description of the Study Area**

The study area (with a population of over 1.7 million and area of 321.03 km2) is the metropolitan area of Tabriz, East-Azerbaijan, Iran, defined by its 10 municipal regions (Figure 1). Tabriz City is located in Azarbaijan geological zone, which is surrounded in the northern region by Alborz, in the south by Semnan and in the west by the Tabriz–Urumiyeh Faults [77]; therefore, it is considered as an area susceptible to earthquakes. Tabriz City continues to the Pontic highlands in Turkey [77,78]. The central and western regions of Iran are comparable with the Azarbaijan geological zone [79], where there are a few important faults [15,78–80]. The Tabriz fault is the most important one near to the city of Tabriz, extending in the northwest–southeast direction from the Zanjan zone and continuing to the northern mountains of Tabriz City [78]. It has been selected because: (1) it is the 5th large city of Iran which is located near the Tabriz fault, and in recent years, the land cover has rarely been updated; (2) it is a hub for cities in the northwest and west of the country (mainly due to better facilities including more job opportunities, higher quality education and more health centers, more people tend to migrate to the city); and (3) for such cities with this level of importance, future land cover prediction and seismic hazard assessment is vital. Seasons in Tabriz are regular, and it has a continental and cold semi-arid climate; at the same time, the average annual precipitation in the study area is around 320 mm, while the average annual temperature is almost 12.6 ◦C [81,82]. The city experiences humid and rainy weather in autumn, while it has a few snowy days during the winter season; at the same time, in spring, the city has a mild and fine climate, and during summers, the region can experience a semi-hot climate [81,82].

**Figure 1.** The geographical extent of the city.

## **3. Materials and Methods**

## *3.1. Database and Data Acquisition*

## 3.1.1. Optical Satellite Data

For the study area, we obtained four cloudless satellite data, from which a Landsat-5 image of the year 1989 was selected as the base image for the study. All the Landsat images were downloaded through the USGS portal. Based on the metadata, only the data with cloud coverage of less than 10% were selected. Then, we searched for and collected RS data between 1989 to 2020, and cloudless Landsat-5 data for 2005 and 2011 were found. Regarding data for the year 2020, an image of Landsat-8 operational land imager (OLI) was collected. The images include seven bands in the range of visible to thermal-infrared for Landsat-5 and nine bands for Landsat-8. The ground resolution of images in optical bands is 30 m, except for a panchromatic band of Landsat-8, which is 15 m. The technical characteristics of the Landsat data are clearly presented in Table 1.

**Table 1.** Characteristics of the optical data used.


## 3.1.2. Synthetic Aperture Radar data

The study was focused on the Landsat missions to carry out the objectives, but the SAR data of Sentinel-1 (SLC and GRD) were also employed for extracting land cover, especially built-up areas, because built-up areas are associated more with seismic hazard assessment; therefore, the validity and reliability of its extraction should be taken into account. Table 2 shows the geometric attributes of SAR data.

**Table 2.** Geometric attributes of SAR data used.


## *3.2. Data Preprocessing and Processing*

First of all, all data of Landsat-5 and 8 were processed for atmospheric, radiometric corrections and the spatial resolution of them was enhanced to 15 m using a panchromatic band of Landsat-8. It is worth mentioning that in pan-sharpening, spectral information will remain unchanged, while the spatial resolution of higher pixel size images will be assigned to the lower one. Here, the 30-m spatial resolution of the Landsat-5 data was enhanced to just 15-m using the panchromatic band of the Landsat-8. These data were imported into TerrSet Software for classification, change detection and prediction using the Markov chain algorithm for the years 2011 and 2030. At the first stage, the regions of interest (ROIs) were extracted carefully; then, by the maximum likelihood algorithm, changes in land cover from Landsat satellite images were detected, classified and mapped. Furthermore, all ancillary data were processed and applied together with classified maps for the prediction steps using the Markov chain model. The overall approach for the current research includes three key procedures: (1) geometric correction of data; (2) classification of optical satellite data and prediction of the future potential landscape of the city; and (3) municipal performance evaluation and seismic hazard assessment (Figure 2). Ancillary reference data were collected from the Municipality of Tabriz

and open street map (OSM) and were applied for training the Markov chain algorithm. These data included a digital elevation model (DEM), buildings, land use, places, railway, roads, green space, waterway and welfare services shapefile of the city (Table 3). Additionally, those data which were collected from the municipality of Tabriz were up to date and were based on the latest changes that occurred in the city.

**Figure 2.** The methodology of the research.


**Table 3.** Information on ancillary reference data.

## *3.3. Model Used in the Study*

#### 3.3.1. Maximum Likelihood

A maximum likelihood classifier was applied to extract surface information from RS data. This defined the statistical values with a normal distribution for each class in image's bands. On the other hand, the algorithm estimated the probability that one pixel would fall into the defined classes. This procedure was continued for all the pixels and the pixels were assigned for those classes that produced the highest probability as follows [83]:

$$\log\_i(\mathbf{x}) - \ln p(\infty\_i) - \frac{1}{2} \ln |\Sigma\_i| - \frac{1}{2} (\mathbf{x} - m\_i) \sum\_i - 1 \ (\mathbf{x} - m\_i) \tag{1}$$

where the number of classes is defined by *gi* (x), which represents the number of imageries' bands, p(∞*i*) describes the probability of class, which occurs in the images, the covariance matrix is defined by *i* ; additionally, the inverse matrix is *<sup>i</sup>* −1, and *mi* is the mean vector [83].

## 3.3.2. Markov Chain

Markov chain is a model from which the future potential landscape can be relatively detected, so that, based on the extracted information from the past data, it detects the future pattern of a land cover [25–27]. In this countable sequence, the chain moves state at discrete time steps [84,85]. It is worth mentioning that this sequence of time process is called a sequence-time Markov chain [84,86]. Markov chains have many applications in different fields. Overall, this model for land cover prediction produces promising findings [4,35]. The model calculated the following formulas [25]:

$$p\_{ij} = \frac{n\_{ij}}{n\_i} \tag{2}$$

$$\sum\_{j=1}^{k} p\_{ij} = 1\tag{3}$$

where transition probability is defined by *pij*, *i* and *j* describe two types of land cover, the total of pixels of each class is shown by *ni* and *nij*, which represents the number of transformed pixels from class *i* to class *j*, and finally, *k* defines the number of land cover classes.

## 3.3.3. Logistic Regression Model

A good deterministic seismic hazard assessment is generally associated with effective and well-approved models [16]. In recent years, the LR model has been widely employed to analyze binary variables and, as a result, it has been introduced as an promising approach in environmental studies [19,42]. Therefore, the LR classification model was adopted in this study. It deals with independent and dependent parameters, where this relationship is nonlinear and can be calculated as Equations (4) and (5) [87,88]:

$$P = \frac{1}{1 + e^{-Z}}\tag{4}$$

where *P* is an earthquake's probability occurrence (0 ≤ *P* ≤ 1), and *e*−*<sup>Z</sup>* is a linear logistic factor (−∞ ≤ *P* ≤ +∞) that is calculated based on Equation (5) [87]:

$$Z = \log \text{ it } (p) = \ln(\frac{p}{1-p}) = b\_0 + b\_1 \mathbf{x}\_1 + \dots + b\_n \mathbf{x}\_n \tag{5}$$

where *Z* is a linear logistic factor, *p* is an earthquake's probability occurrence, *n* is the number of conditioning variables, and *b*<sup>0</sup> is the constant coefficient.

#### Factors Used for Seismic Hazard Assessment

Hazard studies give valuable information about human environments, which, if the results of these studies are taken into account, may protect them from such events in the future. Complex natural hazards such as landslides or flood mapping need considerable data collection and analysis [19,89]. However, in this study, a susceptibility map of potential sites of earthquakes is produced, in which the most important conditioning parameters for it are soil type, proximity to fault lines and lithology condition (Figure 3). A simple probabilistic seismic hazard analysis (PSHA) model was used, which is a useful algorithm, especially when in situ seismic data are not widely available. In the study area, the distribution of the faults is not complex and the Tabriz fault's orientation is straightforward. However, it must be noted that, due to the lack of actual seismic data, this model cannot address uncertainties well [74]. The occurrence probability and intensity of risk assessment depend on selected conditioning factors [17].

**Figure 3.** Conditioning parameters used for susceptibility mapping (**a**) fault lines; (**b**) lithology; and (**c**) soil type.

## *3.4. Accuracy Assessment and Validation*

## 3.4.1. Confusion Matrix for the Classified Maps

Using Google Earth (GE) images, thirty (altogether ninety) ground control points (GCPs) were randomly extracted for each land cover (Figure 4) and then converted into ROIs for accuracy assessment. The overall accuracy and Kappa coefficient were employed in this research. These two models are measured based on Equations (6) and (7), respectively [90,91]:

$$OA = \frac{1}{N} \sum p\_{ii} \tag{6}$$

where *OA* defines the total accuracy of the model, test pixels are described by *N*, and *pii* represents the total number of correctly classified pixels.

**Figure 4.** GCPs for all years.

Kappa coefficient is a statistical model that is employed to measure the reliability of the qualitative items [91]. It is universally accepted that this model is a more robust method than the simple calculations. The Kappa coefficient provides reliable and valuable information for the findings obtained. *OA* must be calculated first in order to measure it.

$$K = (OA - \frac{1}{q})(1 - \frac{1}{q})\tag{7}$$

where *OA* defines the total accuracy of the model; *k* and *q* are Kappa coefficient and unclassified pixels, respectively.

## 3.4.2. Validation of the Predicted Map of the Year 2011

The predicted land cover map using the Markov chain model for 2011 was validated by the generated land cover map of the same year. This was only conducted to determine the reliability and accuracy rate of the model that will be used for the prediction of the future landscape of the year 2030.

## 3.4.3. Validation of Extracted Land Cover Using SAR Data

GRD and SLC products of Sentinel-1 SAR data were used for the validation of mapped land cover. For this reason, a pair of SLC products for two close dates was preprocessed in sentinel application platform (SNAP) software and used for the extraction of land cover using RGB creation in the GIS environment. At the same time, the GRD product was also applied for this matter in order to ensure the complete reliability of the mapped land cover; the reliability of the land cover was essential because it was used to create many maps for the study area.

## **4. Results**

## *4.1. Land Cover Classification*

Reliable detection of landscape change using remote sensing data must strike a balance between affordability and product accuracy [78,92,93]. By applying the maximum likelihood method, vegetation, built-up and bare land surfaces in 1989, 2005, 2011 and 2020 were extracted for the study areas. This information was utilized to create a few maps of buildings by age, municipality performance and seismic hazard assessment. ROI extraction is one the most important steps in land cover classification, from which exact land cover can be extracted, and it also affects the overall accuracy [70,71,94]. Tables 4 and 5 detail the number of pixels and ROI separation characteristics, respectively. The mean pixel count of the extracted ROIs was used for obtaining the spectral signatures of the land cover. Therefore, an image-derived technique was applied for the extraction of the spectral signatures.




**Table 5.** ROI pair separation.

The land cover maps for years 1989, 2005, 2011 and 2020 were generated using the maximum likelihood algorithm. Most of the new built-up areas occurred at the edges of the existing urbanized regions, which are displayed in orange color. Figure 5 details the spatial patterns of classified land cover from 1989 to 2005. The vegetation extent on the map is presented in green, built areas in orange and bare land in light yellow pixels. Figure 6 quantifies the changes which occurred from 1989 to 2005. For the sake of distribution clarity, two column charts of gain/losses and net changes were created for changes that occurred from 1989 to 2005. Moreover, vegetation lost around 20 km2 and gained 18 km<sup>2</sup> from 1989 until 2005 (net change <sup>−</sup>2.63 km2). At the same time, 49.47 km2 was added to the built-up area, while only 4.07 km<sup>2</sup> was removed from it (net change +45.39 km2). Finally, compared to built-up areas, losses for bare land are considerable, so that it lost around 55.35 km2 from its areas and approximately 12 km<sup>2</sup> was added to bare land (net change <sup>−</sup>42.76 km2).

**Figure 5.** Spatial distribution of changes from 1989 to 2005.

**Figure 6.** Statistical attributes for changes which occurred from 1989 to 2005; (**a**) gain and loss (**b**) net changes.

After visual inspection, general information over the city was gathered; the only drawback of our classification using the maximum likelihood algorithm was considering the airport band of Tabriz City as a built-up area. This bias may be because of the similarity of backscatters for roads inside the built-up areas with the airport band; however, this is a negligible area and can be addressed using simple editing using GIS.

Figure 7 represents the spatial trends of land cover from 2011 to 2020, from which increasing vegetation coverage inside the city and beyond is considerable. Two charts of gain/losses and net changes were also provided for changes which occurred from 2011 to 2020 (Figure 8); from these, it can be concluded that from 2011 onwards, only around 4 km2 was added to the built-up areas. For many readers, this should be of great concern, but, based on an interview with the municipality of Tabriz (the interview was performed with a public affairs officer of the municipality on 18/08/2020 through phone call), from 2011 onwards, the city's buildings were constructed and grew vertically, meaning that old buildings with one or two floors were replaced by buildings with more than three floors. A summary of statistical reports for land cover changes from 2011 to 2020 is as follows: (1) net change for vegetation coverage was +20.56 km2, meaning that approximately 33.83 km2 was added to it and 13.27 km2 was subtracted from it; (2) 18.64 km2 was subtracted from built-up areas, while 22.73 km<sup>2</sup> was added to it; and (3) not surprisingly, bare land lost 46.14 km2 and gained 21.68 km2 so that the net change for it can be <sup>−</sup>24.47 km2.

**Figure 7.** Spatial distribution of changes from 2011 to 2020.

Figure 9 clearly shows cross changes from one land cover to the other. For example, those areas that were once built-up areas but were replaced with vegetation coverage are displayed in light green. Dark green represents areas that once were bare land that have become vegetation. The areas indicated by the red color are areas that were originally vegetation but were replaced with built-up areas. At the same time, changes from bare land to built-up areas are indicated by the dark red color. Furthermore, the light yellow color shows areas that were replaced by bare land from vegetation. Additionally, the dark yellow color highlights areas that were built-up areas but then changed to bare grounds. This kind of map is important in showing changes among land cover between two specific years.

**Figure 8.** Statistical attributes for changes which occurred from 2011 to 2020; (**a**) gain and loss (**b**) net changes.

**Figure 9.** Cross change map of land cover.

## *4.2. Future Potential Landscape of Tabriz Using Markov Chain Model*

The land cover potential pattern of the year 2030 was mapped for the study area. To understand the level of reliability of the model used (because there was no information for the year 2030), the land cover map of the year 2011 was also estimated using land cover maps of the years 1989 and 2005. Considering the predicted land cover map of 2011, around 86% of vegetation coverage was forecasted to remain unchanged (which is a high percentage), and changes from built-up area to vegetation were predicted by approximately 1%, while this rate was roughly 4% from bare land to vegetation. It was forecasted that 96% of the built-up areas would remain unchanged, which is also quite high and shows that the model works well, so the prediction for the year 2030 can be reliable to a great extent. The change prediction rates of vegetation to built-up areas and also the bare lands to built-up areas were estimated to be almost 6% and 8%, respectively. Like vegetation, around 86% of bare land was predicted to remain bare land by the year 2011 (which also represents a good prediction rate). Figure 10 and Table 6 detail the spatial pattern and statistical changes in land cover by the year 2011 (using land cover maps of years 1989 and 2005), respectively.

**Table 6.** Probability of land cover changes in 2011 predicted from maps of the years 1989 and 2005.


**Figure 10.** Predicted map of 2011.

Figure 11 and Table 7 highlight the spatial trends and statistical findings of land cover changes by the year 2030, respectively. Considerable findings were extracted regarding the predicted land cover from the map of 2030. Around 74% of vegetation coverage was forecasted to remain unchanged, meaning that almost 26% is likely to be replaced by other types of land cover. Changes from built-up area to vegetation are predicted by approximately 7%, while this rate is quite high for bare land to vegetation, at roughly 26%. It was estimated that 79% of the built-up areas would remain stable as themselves. The change prediction rates of vegetation and bare land to built-up areas are high, at an estimated rate of almost 9% and 11%, respectively. Approximately only 61% of bare land is likely to remain bare land by the year 2030, meaning that based on the changes which occurred until the year of 2020, the municipality plans to change bare land to other types of land cover. One of the considerable results in this regard could be the probability of changes from built-up areas to bare land, at approximately 13% (which is a quite high figure); based on the interview with the municipality, this is maybe because of the reconstruction of buildings over the city that occurred and was recorded by the satellite images used in this study. However, these are only predictions based on changes from the year 2011 to the year 2020.

**Table 7.** Probability of land cover changes in 2030 predicted from maps of years 2011 and 2020.


## *4.3. Building Age Map of the Study Area since 1989*

The built-up areas' distribution by age for the years of 1989, 2005, 2011 and 2020 was extracted (Figure 12 and Table 8). Pixels in light pink color are classified as built-up areas until 1989 and the pink color represents built-up areas that developed after the year of 2005. Areas indicated in red color are defined as newer urban areas that were constructed from 2005 until 2011. The newest built-up areas that have been constructed since 2011 are shown in the dark red color. Most of Tabriz City was constructed before 2011, but relatively new urban areas in and around the study area can be seen, which indicates that urban development has been gradually taking place. The proportion of built-up areas is presented in Table 8, which shows that the total built-up area constructed before 1989 is around 45 km2. In the year 2005, the built-up area doubled to approximately 90 km2. Almost 21 km2 was added to built-up areas by the year 2011. Not surprisingly, only approximately 4 km2 has been added to the built-up areas since 2011 (this does not mean that urbanization has stopped since then); this is because buildings have been reconstructed vertically (a few floors) instead of containing only one or two floors. Our interview with the municipality also confirmed that the old buildings with one or two floors are being reconstructed and replaced with buildings with three or more floors. This has

good advantages, such as providing more land with the municipality for establishing other projects, while there is sufficient housing for citizens as well.

**Figure 12.** Spatial distribution of building age.

**Table 8.** Area of constructed regions by four different years.


Urbanization Rate

Based on the built-up area extracted from the classified maps in this study, urbanization rate (UR) was calculated using the following user-defined equation:

$$\text{LIR} = \frac{A}{T} \tag{8}$$

where *UR* is the urbanization rate, *A* is an extended area of built-up areas for each period, and *T* contributes to the time passed for urban growth. According to the equation, *UR* for a different period related to the current study is as follows: it is worth mentioning that the results are km2 per year.

$$\text{LIR}(2005) = \frac{45.68}{16} = 2.85$$

$$\text{LIR}(2011) = \frac{21.72}{6} = 3.62$$

$$UR(2020) = \frac{4.11}{9} = 0.45$$

## *4.4. Municipal Performance Evaluation for 10 Municipal Zones of Tabriz City*

Based on the changes that occurred from the year 2011 to 2020, the performance of municipal zones of Tabriz city was evaluated (considering the changes towards more built-up areas and green space along with less bare land), meaning that when more bare land for a municipal zone was converted to built-up and vegetation areas, it was considered that the zone worked well. Following our evaluation, the best and worst municipalities are municipal numbers one and four, respectively. Figure 13 and Table 9 show spatial and statistical changes in different types of land cover for each municipality, respectively.


**Table 9.** Quantitative results for municipal performance evaluation since 2011.

**Figure 13.** Municipal performance evaluation map.

## *4.5. Seismic Hazard Assessment*

Four susceptible zones were finally reclassified. Figure 14 clearly shows the spatial patterns of areas to susceptible to earthquakes concerning the municipal zones of the city. Most of the municipal zone numbers 1, 9 and 10 are located in the low susceptibility zone, while the entire municipal zone number 8 and most of the municipal zone numbers 3, 4, 6 and 7 are located in the very high susceptibility zone.

**Figure 14.** Earthquake susceptibility mapping.

## *4.6. Accuracy Assessment and Validation for This Study*

## 4.6.1. Confusion Matrix for the Classified Maps

The land cover classification mapping using RS data is a relatively easy effort, while the accuracy and inaccuracy of it depends on proper data and models used [2,61]. An accuracy assessment using the confusion matrix was accomplished for all four classification maps. The overall accuracy for the classified maps of 1989, 2005, 2011 and 2020 was around 96% 96%, 93% and 94%, respectively. Meanwhile, these values for the Kappa coefficient were almost 93%, 92%, 85% and 88%, respectively (Table 10).


**Table 10.** Statistical results using the confusion matrix for all classified maps.

## 4.6.2. Validation of the Predicted Map of the Year 2011

Using the land cover map of 2011, the predicted map of the same year was validated. Since there was not any information for the year 2030, the validation was not possible. However, fortunately, the prediction map of the year 2011 was well validated (meaning that it only predicted a few areas incorrectly; most were predicted well); therefore, it can be concluded that the prediction map of the year 2030 can be also correct to a great extent. These land cover products were then used to validate urban extent extraction, which confirmed that land cover extraction was done successfully. Validation interpretation was based on these three data: (1) initial land cover was the classified map of the year 2005 (2); predicted land cover for the year 2011; and (3) validation land cover (classified map of the year 2011). Figure 15 represents the validation of the predicted map of the year 2011. However, for interpretation of the image, two examples are presented here: (1) 1/1/1 means that these areas in all aforementioned images were vegetation, and (2) 2/3/3 means that these areas were bare land originally but were predicted as built-up areas.

4.6.3. Validation of the Extracted Land Cover for the Year 2020 as a Basis for Seismic Hazard Assessment

The magnitude of errors using conventional methods is a complex issue from which the extracted land cover from them cannot be directly applied for understanding the changes which occurred [3,52]. Additionally, uncertainties are inherent aspects of remotely sensed studies [3]; to minimize uncertainties, SAR data were also applied. This attempt ensured that urban land cover was not missed using optical images. To successfully validate the extracted built-up areas which were utilized for seismic hazard assessment of the city, SAR data were also used. A few small areas were marked by a few geometric symbols (square-shaped) in each set of satellite data (Figure 16). These areas were then enlarged and displayed with different shapes for each type of land cover, comparing the built-up surface from two SAR data of GRD and SLC products, ensuring that the built-up areas in both were successfully matched. After preprocessing and processing of SAR data in the SNAP environment, both bands of VV and VH were employed for RGB creation in the GIS environment. The VV band of the slave imagery was used for R and the VH band of the master one was employed for G and B windows.

**Figure 15.** Validation of the predicted map of the year 2011.

**Figure 16.** Land cover maps using both optical and SAR data.

#### **5. Discussion**

Generally, a few diverse factors may affect the results of change detection and prediction [9,84]. The land cover's spatiotemporal pattern and characteristics in this study were completely different from those which have been measured before, in which a comprehensive study including land cover change detection, future potential landscape, distribution of buildings by age, municipal performance evaluation and building damage assessment was carried out for the metropolitan area of Tabriz, especially for the municipality of Tabriz that has suffered from a lack of such studies. Previous studies were mainly focused on simple methods for measuring the land cover, seismic hazard and building vulnerability for the study area [80,95,96]. During the first period (1989–2005), the built-up area increased as a result of population growth and migration to the city, which was mainly based on destroying vegetation and bare land. During the second period (2011–2020), bare land was replaced by other types of land cover and, in this period of time, one of the most considerable findings was increased vegetation across the city, which, based on the findings, reflects the efforts of the municipal regions to increase the vegetation to a satisfactory level. Urban growth was mainly observed in the bare land and the vegetated areas, far from the economical areas.

The dynamics of the land cover are correlated [77]. Here, built-up areas grew considerably in the first period (1989 to 2005), while in the second period (2011 to 2020), this growth was negligible. For the second period, we found that vegetation areas experienced more positive changes than the first period among land cover. The built-up areas in both periods (1989 to 2005 and 2011 to 2020) showed the largest degree of change among all land covers, which could be linked to the rapid urbanization in Tabriz. The change in speed of bare land was relatively fast in both periods, in which it underwent the most significant land cover changes for the entire period. In addition, vegetation was directly linked to the civil projects of the municipality that turned bare land and old built-up areas into green lands. This analysis suggests that the land cover in Tabriz has considerably changed during the last three decades.

According to the findings for the year 2030, the general trend of change is toward more vegetation and built-up areas as well as less bare land, meaning that the municipality plans to convert more bare land to other types of land cover. Specifically, in the year of 2030, the vegetation and built-up areas will preserve most of their areas and will be larger because of the change from bare land to these types of land cover. Besides this, almost 60 percent of the area of bare land will be preserved, because always bare land will be used for the development of new projects. This implies that the landscape pattern of Tabriz has a tendency to be more optimized in this period.

Based on the results of the urbanization rate, from the year 1989 until 2005, the city of Tabriz experienced growth in urban areas of around 2.85 km<sup>2</sup> per year. At the same time, this rate was 3.62 km2 per year from 2005 until 2011. Additionally, from 2011 to the year 2020, this rate was only 0.45 km<sup>2</sup> annually.

Seismic hazard assessment is always an essential part of sustainable development projects for urban areas [7,15]. Seismic hazard assessment produces a greater cost efficiency when focusing only on urban areas (rather than the entire area including vegetation and bare land), which has the greatest impact on people when they are destroyed. The application of seismic hazard assessment has not yet been conducted for a city like Tabriz, which is growing fast based on the population rate and its situation as a hub to the other cities of the region. Therefore, a seismic hazard assessment (even a simple one) could enable the local authorities and the policy-makers to direct urbanization to those areas with low or moderate susceptibility. Concerning the earthquake risk for the city, policy-makers should take this into account for future urban sprawl, meaning that they can design stricter policies for new buildings that are constructed or reconstructed in more susceptible areas.

Model validation is important to assess the level of the models' reliability and validity. Validation of the created maps was performed at three stages. In the first, the classified outputs of the years 1989, 2005, 2011 and 2020 were validated using a confusion matrix. In the second step, and since the land cover map of the year 2020 was selected for damage assessment, two other maps from Sentinel-1 SAR

data were also employed. Finally, a classified map of the year 2011 was used to validate the prediction map of the year 2011 (only to assess whether the model could predict the future potential landscape or not, so that its prediction for the year 2020 could also be considered correct).

Based on the findings of this study, for arid and semi-arid regions (like the current study area), the maintenance of the existing vegetated areas rather than planting more grasses and trees in less suitable areas is recommended. Although this study, as well as many previous works, has demonstrated that the remotely sensed data and techniques can be well applied for monitoring changes in cities, we recommend that more high-resolution satellite imageries be used to gain further insights into such changes. Future research should focus on the deep learning techniques for change detection and the prediction of land cover; more details on seismic and risk assessments can also be obtained using deep learning algorithms for the study area. The only major limitation of this study was encountered when obtaining ancillary data from the municipality of Tabriz.

## **6. Conclusions**

Cities need comprehensive and innovative plans in order to ensure progress based on sustainable development. Although it is very difficult to obtain absolute results from remotely sensed data, relative findings can be captured, which can be effective for any future planning. This study has emphasized changes in land cover and the future landscape in Tabriz City. Other important issues that the current research was focused on include information on building age, municipal performance evaluation and building damage assessment, which contributes to earthquake damage estimation. This study has also compared the results of optical satellite imagery with SAR data to extract the spatial distribution of buildings for the year 2020, which was the base map to evaluate municipal zone performances and seismic hazard assessment. The main findings of the current study are as follows: Landsat images for the years 1989, 2005, 2011 and 2020 were used to quantify the land cover changes from 1989 until 2020 and the results using the confusion matrix were promising. At the same time, by using and comparing SAR data, the accuracy of built-up areas for the year 2020 was well validated and verified. Referring to the assessment of the distribution of built-up areas by age for Tabriz City, we found that most of the built-up areas had been developed before 2011, and from then onwards, the city has been progressing vertically. Seismic hazard assessment for the future of the city was conducted by using a logistic regression model, from which results indicated that municipal zones 1 and 9 are located inside low susceptibility areas, while municipal zones 4, 6, 7, 8, and also most of zones 3 and 10, are located in highly susceptible regions. Further findings revealed that land cover prediction by using the Markov chain model provided a good opportunity to identify the future potential landscape of the city. Finally, based on the land cover maps of 2011 and 2020, the performances of the municipal zones were evaluated, from which results showed that municipal zone 1 followed by zone 5 have the best performances among all. Besides this, the performance of municipal zone 4 is negligible, as is much of municipal zone 6.

**Author Contributions:** A.M., S.K., K.V.K. and M.M. collected data and finalized the study; A.M., S.K. and M.M. worked on the methodology; A.M. and S.K. performed software analysis and validated the research findings; A.M. and S.K. finalized the formal analysis; A.M., S.K. and M.M. wrote the original draft; S.K., K.V.K. and M.M. reviewed and edited the original version; S.K., K.V.K. and M.M. funded the research. All authors have read and agreed to the published version of the manuscript.

**Funding:** The presented work is supported by National Elites Foundation of Iran (No. 102/1670) and the Japan Society for the Promotion of Science (JSPS), KAKENHI number 20H02411.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
